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Preface 
Preface to "A Brief Introduction to Engineering Computation with 
MATLAB" 


In my tenth year at the Institute, I dedicate this book to the BCIT 
community. 


The primary purpose of writing a book and distributing it free-of-charge is 
to extend my gratitude to BCIT. I am particularly thrilled to do it with this 
textbook because it is a product of many learning opportunities BCIT has 
offered me over a period of several years. What follows is a brief 
background on how this book came to be. 


My post-secondary teaching career began on 22 January 2001 at the Pacific 
Marine Training Campus of BCIT when I logged on to a Unix workstation 
to instruct in the Propulsion Plant Simulator. That has been a major 
milestone in many ways in my professional life. While learning inner 
workings of Unix operating system (OS), I also made a discovery and that 
discovery profoundly changed my view on how I thought the world 
operated. The discovery was the GNU/Linux OS and open source software 
(OSS) movement through several books, most notably Just for Fun: The 
Story of an Accidental Revolutionary[footnote] and The Cathedral and the 
Bazaar[footnote]. I was convinced that the collective power of connected 
individuals around the world and the global infrastructure of the Internet 
had the potential to change the ways the world functioned. 

Just for Fun: The Story of an Accidental Revolutionary by L. Torvalds and 
D. Diamond, New York: HarperCollins Publishers. © 2001 

The Cathedral and the Bazaar by E. S. Raymond, Sebastopol: O’Reilly 
Media. © 1999 


In the last 10 years, BCIT has allowed me to study various subjects through 
its Professional Development (PD) programs for which I am very grateful. I 
learned a great deal in PD courses and in one of the recent ones, I had two 
déja vu moments similar to my discovery of OSS movement. The first one 
occurred when I began reading The Wealth of Networks|[footnote] and the 
second one when I found about Connexions. The former was a confirmation 
of my 10-year old discovery and the latter is what I am using to write this 
book. Connexions is a web-based curricular content authoring and 


publishing technology that I believe has a growing potential for writing and 
distributing free-of-charge learning materials. 

The Wealth of Networks by Y. Benkler, New Haven: Yale University Press. © 
2006 


Thus, motivation for this book stems from the notions that were generated 
by the OSS movement. The book was written to pay a small token of 
appreciation to BCIT and I hope it will be a contribution to the open 
educational resources repository. 


Serhat Beyenir 
North Vancouver, B. C. 
25 October 2011 


Study Guide 
Study guide for learners. 


MATLAB, a sub-course of Computer Technology 1 and this text are 
specifically designed for students with no programming experience. 
However, students are expected to be proficient in First Year Mathematics 
and Sciences and access to good reference books are highly recommended. 
I also assume that students have a working knowledge of the Mac OS X or 
Microsoft Windows operating systems. 


The strategic goal of the course and book is to provide learners with an 
appreciation for the role computation plays in solving engineering 
problems. The MATLAB specific skills that I would like students to acquire 
are as follows: 


e Write scripts to solve engineering problems including interpolation, 
numerical integration and regression analysis, 

e Plot graphs to visualize, analyze and present numerical data, 

e Publish reports. 


The best way to learn about engineering computation is to actually do it. We 
will therefore solve many engineering problems mainly using a recent 
version of MATLAB in this book. Since the primary focus is engineering 
computation, we will concentrate on the mathematical solutions and, to a 
limited extent, the graphical user interface (GUI) features of MATLAB. 


Learning a new skill, especially a computer program, can be an 
overwhelming experience. To make the best of this process, students are 
encouraged to observe the following guidelines that have proven to work 
well: 


e Plan to study 2 hours outside of class for every hour inside of class, 

e Practice, practice, practice: As the old saying goes, practice makes one 
perfect or perhaps we should modify that statement: Good practice 
makes one perfect, 

e Buddy system: Study with a classmate. Helping one another drastically 
improves your understanding of the material. Particularly, students are 
advised to work the problem sets in this fashion, 


e Muddy points: Make a note of muddy points as they may occur during 
lectures and email your notes to me. I will address those issues at the 
beginning of the next class, 

e Open book exam: Do not try to memorize commands, functions or 
their syntax but learn where and how to find that information. Through 
many exercises and problem sets you will have solved by the end of 
the course, most computational routines will become second nature to 
you. The exam is open book, so keep your learning materials and m- 
files well organized. 


What is MATLAB? 
A brief introduction to MATLAB. 


Why MATLAB 
GUI 
he Ip 


Useful Commands 
jp and Functions 
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MATLAB stands for MATrix LABoratory (see wikipedia) and is a 
commercial software application written by The MathWorks, Inc. When 
you first use MATLAB, you can think of it as a glorified calculator 
allowing you to perform engineering calculations and plot data. However, 
MATLAB is more than an advanced scientific calculator, for example 
MATLAB's sophisticated numerical computation environment also allows 
us to analyze data, simulate engineering systems, document and share our 
code with others. 


Why Use MATLAB? 


MATLAB has become a defacto standard in many fields of engineering and 
science. Even a casual exploration of MATLAB should unveil its 
computational power however a closer look at MATLAB's graphics and 
data analysis tools as well as interaction with other applications and 
programing languages prove why MATLAB is a very strong application for 
technical computing. 


The standard MATLAB installation includes graphics features to visualize 
engineering and scientific data in 2-D and 3-D plots. We can interactivity 
build graphs and generate MATLAB command output that can be saved for 
use in the future. The saved-instructions can be called again with different 
data set to build new plots. The plots created with MATLAB can be 
exported in various file formats (e.g. .jpg, .png) to embed in Microsoft 
Word documents or PowerPoint slideshows. 


MATLAB also contains interactive tools to explore and analyze data. For 
example, we can visualize data with one of the many plotting routines, 
zoom in to plots to take measurements, perform statistical calculations, fit 
curves to data and evaluate the obtained expression for a desired value. 


MATLAB interacts with other applications (e.g. Microsoft Excel) and can 
be called from C code, C++ or Fortran programming language. 


Running MATLAB 


To use MATLAB, it must be installed on your computer and you can start it 
just like you start any application on your system or you must have access 
to a network where it is available. 


BCIT holds a Total Student Headcount (TSH) license for Mathworks 
software and this allows students to install MathWorks software on their 
personally-owned computers. 


Matlab download and installation instructions for BCIT students 


For your install, choose the latest version available and install MATLAB, 
Curve Fitting Toolbox and Symbolic Math Toolbox. 


In addition, MATLAB Online provides access to MATLAB from your web 
browser. Just log in to use MATLAB: 


Access MATLAB From Your Web Browser 


The MATLAB Desktop 


When you start the MATLAB program, it displays the MATLAB desktop. 
The desktop is a set of tools (graphical user interfaces or GUIs) for 
managing files, variables, and applications associated with MATLAB. The 
first time you start MATLAB, the desktop appears with the default layout, 
as shown in the following illustration. 


Lei Eka 
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Command Window 
fe >> cid| 
J) SRECYCLE.BIN 
J) Temporaryltems 
J» 2008 
builds 
) Desktop 
4 Documents 
LabVIEW Data 
MATLAB 
.) MyCourses 
}) PersonalFiles 
» public_html! CommandHistoy ss © 
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Select a file to view details 


The MATLAB Desktop. 


Command Window 


The Command Window is where we execute MATLAB commands. We 
enter statements at the Command Window prompt. The prompt can be any 


one of the following: 


e Trial>> indicates that the Command Window is in normal mode and 
the MATLAB license will expire after the trial period ends. 

e EDU>> indicates that the Command Window is in normal mode, in 
MATLAB Student Version. 

e >> indicates that the Command Window is in normal mode. 


Command Window 


This is a Classroom License for instructional use only. 
Research and commercial use is prohibited. 


fe >> 


The Command Window. 


Command History 


The Command History is a log of the commands we have executed in the 
command window. 


Command History 
%--— 03/10/2014 12:02... 


The Command History. 


Workspace 


The workspace consists of a set of variables stored in memory during a 
MATLAB session. To open the Workspace browser, select Desktop > 
Workspace in the MATLAB desktop, or type 


>> workspace 


at the Command Window prompt. 


Workspace. 


Current Folder 


The Current Folder is like the Finder in Mac OS X or Windows Explorer in 
Windows operating systems and allows us to browse through the files and 
folders. The Current Folder also displays details about files in your current 
directory and within the hierarchy of the folders it contains. 


@o> 8 Wau: +p 


Current Folder. 


. Documents 
)) LabVIEW Data 


) PersonalFiles 
public_html 
» slprj 
J) temp 
) untitled_grt_rtw 
J) UserProfile 
J userprofiler2 


Current Folder docked on 
the desktop. 


Tool Strip 


The tool strip contains global tabs, Home, Plots and Apps. Contextual tabs 
become available when you need them. 


[aoe 


Tool Strip. 


Tabs. 


The plots tab allows us to plot various types of graphs quickly and easily. 


Q)| Search Documentation SiS 
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The plots tab. 


The apps tab gives quick access to interactive applications within MATLAB 
environment. 


The apps tab. 


Layout button allows us to change the desktop layout or go back to the 
default configuration. 


Layout 
button. 


Toolbar 


The MATLAB toolbar provides on-screen buttons to access frequently used 
features such as, copy, paste, undo and redo. 


Toolbar. 


Keyboard shortcuts 


MATLAB provides keyboard shortcuts for viewing a history of commands 
and listing contextual help. 


1. The up arrow key, 
2. The tab key, 
3. The semicolon symbol. 


The Up Arrow Key 

Suppose we want to enter the following equation: 
>> y=sin(45) 

But we mistakenly entered 

>> y=sine(45) 

MATLAB returns the following prompt: 


22? Undefined function or method ‘'sine' for input 
arguments of type 'double'. 


Instead of retyping the equation, press the up arrow key, the mistakenly 
entered line is displayed. Using the left arrow key, move the cursor to the 
misspelled letter. Make the correction and press Return or Enter to execute 
the command. 


Pressing the up arrow key repeatedly recalls the previously entered 
commands. Likewise, typing the first characters of previously entered line 
and pressing the up arrow key displays the full command line. To execute 
that line, simply press the Return or Enter key. 


The Tab Key 


Suppose you forgot how to enter the square root command. Begin typing 
Y=sq in the command prompt: 


>> y=sq 


Then press the tab key and scroll down to sqrt. Select it and press Return 
or Enter key. 


>> y=sqrt 


The Semicolon Symbol 


The semicolon symbol at the end of a line suppresses the screen output. 
This is useful when you want to keep your command window clean. 


Type the following entry and press the Return key: 
>> y=2+2 
The following output is displayed: 
v= 
4 
Now, press the up arrow key to recall our initial entry 
>> y=2+2 
And insert a semicolon as follows: 
>> y=24+2; 
No numerical result is displayed however MATLAB stores the value of y in 


the memory. We can recall the value y by simply typing y and pressing 
Return. 


MATLAB Help 


MATLAB comes with three forms of online help: help, doc and demos. 


Help 


Typing help in the Command Window lists all primary help topics. You can 
display a topic by clicking on the link. 


>> help 


ERE 


> H: > MATLAB > 


|/2 4. Chp2_Exercise9 
|=  Chp3_Exercisea 
\]2 L. EngineeringComputati... 
J figs 
& } html 
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No details available 


Bi 


>> help 
HELP topics: 


Nn) 6.75 


matlab\testframework 


matlab\demos 
Matlab\graph2d 
matlab\graph3d 
Matlab\graphics 
matilab\plottools 
Matlab\scribe 
Matlab\specgraph 
Matlab\uitools 
toolbox\local 
matlab\optimfun 
Matlab\codetools 
matlab\datafun 
Matlab\datamanager 
Matlab\datatypes 
matlab\elfun 
matlab\elmat 
matilab\funfun 
matlab\general 
matlab\gquide 
matilab\helptools 
matlab\iofun 


marlah\lanc 


(No table of contents file) 

(No table of contents file) 
Examples. 

Two dimensional graphs. 

Three dimensional graphs. 

Handle Graphics. 

Graphical plot editing tools 
Annotation and Plot Editing. 
Specialized graphs. 

Graphical user interface compone 
General preferences and configuz 
Optimization and root finding. 
Commands for creating and debuge 
Data analysis and Fourier tranat 
(No table of contents file) 

Data types and structures. 
Elementary math functions. 
Elementary matrices and matrix n 
Function functions and ODE solve 
General purpose commands. 
Graphical user interface design 
Help commands. 

File input and output. 


Procramming landace constructs 


[1,2,3,4,5,6,7,8, 
[0.5403,-0.4161,- 


Command History 
‘cic 


help 


help 
cle 
help clic 
help 
clc 
help 


Help. 


Or if you know the command or function you need help with, you can type 
help followed by the command or function. For example to learn about c1c 
command, type help clc at the command prompt: 


>> help clc 


Patan Ea 
0 ae - 


>> help cic 
elc Clear command window. r 11,2,3,4,5,6,7.8 
ele clears the command window and homes the cursor. 
[0.5403,-0.4161,- 


See also home. 


Reference page in Help browser 


doc cl 
_) M-Files_HeatTransfer mene 


Db siprj 

) _.DS_Store 
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©) AcetyleneBottle.m 

1) AcetyleneBottlelnterac... 

3) AcetyleneBottlelnterac... 

Ub Cho? Fxercise9.7in | y=cos (x) 
Chp3_Exercise4 (Folder) plot (y) 


help 


No details available 


help 
clic 


help clic 


The output of >> help clc command. 


Also try the following command: >> help clear 


: 


Se i > Hu: > MAT 


>> help clear 
2B Chp2_Exercise9 clear Clear variables eis functions from memory. r] [1,2,3,4,5,6,7,8,9,, 
clear removes all variables from the workspace. 
a mY Chp3_Exercise4 clear VARIABLES does the same thing. (0.5403,-0.4161,- 
» EngineeringComputati... clear GLOBAL removes all global variables. 
& figs clear FUNCTIONS removes all compiled MATLAB and MEX-functions 


J htm 
) M-Files_HeatTransfer 
» slprj 
() DS Store clear IMPORT clears the base import list. It can only be iss 
ia DS Store command prompt. It cannot be used in a function. 


meee clear CLASSES is the same as clear ALL except that class defi 
© AcetyleneBottlelnterac... are also cleared. If any objects exist outside the workspace Command History 
© AcetyleneBottlelnterac... userdata or persistent in a locked program file) a warning wi Ee 
1) Chn? FExercise9.7in issued and the class definition will not be cleared. clear Ci 
Chp3_Exercised (Folder) be used if the number or names of fields in a class are chanc 


clear ALL removes all variables, globals, functions and MEX 1] 
clear ALL at the command prompt also clears the base import 1 


doc 
demos 
clc 
clear JAVA is the same as clear ALL except that java classes help 
dynamic java path (defined using JAVACLASSPATH) are also cle cle 
help cic 
clear VAR1 VAR2 ... clears the variables specified. The wilde help 
character '*"' can be used to clear variables that match a pat 
instance, clear X* clears all the variables in the current we 
that start with X. 


No details available 
cle 


help 
clc 


help clear 


The output of >> help clear command. 


To learn about sine function, type help sin at the command prompt: 


>> help sin 


Doc 


Obviously, to use he 1p effectively, you need to know what you are looking 
for. Often times, especially when you first start learning an application, it is 
usually difficult to ask the right questions. In the case of MATLAB, doc 
command is generally better than he Lp. If you type doc in the command 
prompt, MATLAB opens a browser from where you can obtain help easier: 


>> doc 


@ Help 


@% 3 - G7 | Home | + 


Search Documentation 


Installation Release Notes 


MATLAB 

Simulink 

Control System Toolbox 
Data Acquisition Toolbox 
DSP System Toolbox 
Embedded Coder 

Image Processing Toolbox 
Instrument Control Toolbox 
MATLAB Coder 

MATLAB Report Generator 
Optimization Toolbox 
Real-Time Windows Target 
Signal Processing Toolbox 


SimDriveline 


SimElectronics 


SimEvents 

SimHydraulics 
SimMechanics 
SimPowerSystems 
Simscape 

Simulink 3D Animation 
Simulink Coder 

Simulink Control Design 
Simulink Design Optimization 
Simulink Real-Time 

Simulink Report Generator 
Stateflow 

Symbolic Math Toolbox 
System Identification Toolbox 


Built-in MATLAB Documentation. 


Like using help sin, try typing doc sin in the command prompt: 


>> doc sin 


Demos 


You can learn more about MATLAB through demos by typing demo in the 
command prompt, a list of links to demos will open in Help Browser. 
Demos and online seminars are available at product demos and online 


seminars. 


>> demo 


» 3 xk- MATLAB Examples x| +] Bo elo) ~ 
Contents z Search Documentation pe] 
Documentation Center e MATLAB 
v MATLAB 


> Getting Started with MATLAB 
MATLAB Examples 
Release Notes : 
Functions On this page... 


> Language Fundamentals 
> Mathematics 


> Graphics | Language Fundamentals 
> Programming Scripts and Functions | Mathematics 

> Data and File Management | Graphics 

> GUI Building 


Programming 
> Advanced Software Development eae 
| Data Import and Expo 


> Desktop Environment 


Creating Graphical User Interfaces 
| Other Examples 


| New Features Videos 


Getting Started More Examples 


em Getting Started with MATLAB (5 min, 6 sec) EB Video 


SS Working in The Development Environment (5 min, 21 sec) EB Video 


Language Fundamentals 


Create a Structure Array @ Script 


Built-in MATLAB Demos. 


Useful Commands and Functions 


For a detailed explanation and examples for each of the following type 
‘help function’ (without quotes) at the MATLAB prompt. 


Command/Function 
cle 

clear 

who, whos 
workspace 

cd 


pwd 


computer 


ver 


quit 


exit 


Summary of Key Points 


Meaning 

Clear Command Window 
Remove items from workspace 
List variables in workspace 
Display Workspace browser 
Change working directory 
Display current directory 


Identify information about computer on 
which MATLAB is running 


Display version information for 
MathWorks products 


Terminate MATLAB 


Terminate MATLAB (same as quit) 


Useful commands and functions 


1. MATLAB is a popular technical computing application and 
MathWorks offers a trial version of MATLAB on their website, 

2. The MATLAB Desktop consists of Command Window, Command 
History, Workspace, Current Folder and Start Button, 

3. The up/down arrow keys, the tab key and the semicolon are convenient 
tools to use the Command Window, 

4. MATLAB features an online help, doc and demo, 


5. Various commands and functions make MATLAB experience easier, 
forexample,clc, clear and exit. 


Problem Set 
Problem Set for What is MATLAB? 
Exercise: 


Problem: Learn about the following terms using he 1p command: 


1. workspace 
2. plot 

3. clear 

4. format 

9. roots 


Solution: 


1.>> help workspace WORKSPACE Open Workspace 
browser to manage workspace WORKSPACE Opens 
the Workspace browser with a view of the 
variables in the current Workspace. 
Displayed variables may be viewed, 
manipulated, saved, and cleared. See also 
whos, openvar, save. Reference page in Help 
browser doc workspace >> 

2.22 help plot PLOT Linear plot. PLOT(X,Y) 
plots vector Y versus vector X. If X or Y is 
a matrix, then the vector is plotted versus 
the rows or columns of the matrix, whichever 
Tine Ups ii x, 1s) a) scalar and ¥ 1S a Vector, 
disconnected line objects are created and 
plotted as discrete points vertically at X. 

3.>> help clear CLEAR Clear variables and 
functions from memory. CLEAR removes all 
variables from the workspace. CLEAR 
VARIABLES does the same thing. CLEAR GLOBAL 
removes all global variables. CLEAR 


FUNCTIONS removes all compiled M- and MEX- 
functions. CLEAR ALL removes all variables, 
globals, functions and MEX links. CLEAR ALL 
at the command prompt also removes the Java 
packages AMPOre PiSt tava. toon 

.>> help format FORMAT Set output format. 
FORMAT with no inputs sets the output format 
to the default appropriate for the class of 
the variable. For float variables, the 
default aS FORMAT SHORT: .: 2... 

.>> help roots ROOTS Find polynomial roots. 
ROOTS(C) computes the roots of the 
polynomial whose coefficients are the 
elements of the vector C. If C has N+1 
components, the polynomial is C(1)*X‘N + 

fe (CN) 2 Xe ONG Ih) Sees care 


Essentials 
Essential skills to use MATLAB effectively. 


Basic Computation 
Functions 


Variables 
p Linear Equations 


Pol yromals 
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Learning a new skill, especially a computer program in this case, can be 
overwhelming. However, if we build on what we already know, the process 
can be handled rather effectively. In the preceding chapter we learned about 
MATLAB Graphical User Interface (GUI) and how to get help. Knowing 
the GUI, we will use basic math skills in MATLAB to solve linear 
equations and find roots of polynomials in this chapter. 


Basic Computation 


Mathematical Operators 


The evaluation of expressions is accomplished with arithmetic operators as 
we use them in scientific calculators. Note the addtional operators shown in 
the table below: 


Operator 


+ 


Operators 


Name 


Plus 


Minus 


Asterisk 


Forward Slash 


Back Slash 


Caret 


Dot Asterisk 


Dot Slash 


Dot Back 
Slash 


Dot Caret 


Description 
Addition 
Subtraction 
Multiplication 
Division 

Left Matrix Division 
Power 


Array multiplication (element- 
wise) 


Right array divide (element-wise) 


Left array divide (element-wise) 


Array power (element-wise) 


Note:The backslash operator is used to solve linear systems of equations, 


see [link]. 


Note: Matrix is a rectangular array of numbers and formed by rows and 


i 2 et 
5 6 7 8 
columns. For example A = . In this example A 
9 10 11 12 
13 14 15 16 


consists of 4 rows and 4 columns and therefore is a 4x4 matrix. (see 
Wikipedia). 


Note:Row vector is a special matrix that contains only one row. In other 
words, a row vector is a 1xn matrix where n is the number of elements in 
the row vector. B= (1 2 3 4 5) 


Note:Column vector is also a special matrix. As the term implies, it 
contains only one column. A column vector is an nx1 matrix where n is the 
1 


number of elements in the column vector. C = 


ao ® Ww bh 


Note:Array operations refer to element-wise calculations on the arrays, for 
example if x is an a by b matrix and y is ac by d matrix then x.*y can be 

performed only if a=c and b=d. Consider the following example, x consists 
of 2 rows and 3 columns and therefore it is a 2x3 matrix. Likewise, y has 2 


12-3 
rows and 3 columns and an array operation is possible. x = ( ) 


4 5 6 

10 20 30 10 40 90 

and y = then x.*y = 
40 50 60 160 250 360 


Example: 


The following figure illustrates a typical calculation in the Command 


Window. 


Current Folder 


} SRECYCLE.BIN 
) Temporaryltems 
» 2008 
D builds 

.) Custom Office Templates 
) Desktop 
). Documents 
) LabVIEW Data 
MATLAB 
) MyCourses 

). PersonalFiles 


Select a file to view details 


a’ a 
SP? EAs H> 


Command Window 
>> 242 


>> 40/5 


ans = 


>>-2"3 


Operator Precedence 


|_| 

Command History © 

x=12 4| 

y=55 

z=x+y 

buried_duct 

Problem2014 2 

1 
2x 10 

cle 

2+2 

5*8 

40/5 

2*3 


%-- 03/06/2015 09:20 --% 


Basic arithmetic in the command window. 


MATLAB allows us to build mathematical expressions with any 
combination of arithmetic operators. The order of operations are set by 
precedence levels in which MATLAB evaluates an expression from left to 
right. The precedence rules for MATLAB operators are shown in the list 


below from the highest precedence level to the lowest. 


1. Parentheses () 

2. Power (‘) 

3. Multiplication (*), right division (/), left division (\) 
4. Addition (+), subtraction (-) 


Mathematical Functions 


MATLAB has all of the usual mathematical functions found on a scientific 
calculator including square root, logarithm, and sine. 


Note: Typing 

pi 

returns the number 3.1416. To find the sine of pi, type in 
sin(p1) 


and press enter. 


Note:The arguments in trigonometric functions are in radians. Multiply 
degrees by pi/180 to get radians. For example, to calculate sin(90), type in 


sin(90*pi/180) 


Note:In MATLAB 


log 


returns the natural logarithm of the value. To find the In of 10, type in 
log(10) and press enter, (ans = 2.3026). 


Note:MATLAB accepts 
log10 


for common (base 10) logarithm. To find the log of 10, type in log10(10) 
and press enter, (ans = 1). 


Practice the following examples to familiarize yourself with the common 
mathematical functions. Be sure to read the relevant he 1p and doc pages 
for functions that are not self explanatory. 


Example: 
Calculate the following quantities: 
93 
iP 32-1? 
50 
oF =? for d=2 


MATLAB inputs and outputs are as follows: 


aks —— is entered by typing 243/( 342-1) (ans = 1) 
2. 59-5 — J is entered by typing sqrt (5) -1 (ans = 1.2361) 


3. <a for d=2 is entered by typing pi/4* 2/2 (ans = 3.1416) 


Example: 
Calculate the following exponential and logarithmic quantities: 


ie 
2; In(5"°) 
3. log 10° 


MATLAB inputs and outputs are as follows: 


1.exp(2) (ans = 7.3891) 
2. log( (5410) ) (ans = 16.0944) 
3. 1log10(1045) (ans =5) 


Example: 
Calculate the following trigonometric quantities: 


1. cos(Z) 
2. tan(45) 
3. sin(7) + cos(45) 


MATLAB inputs and outputs are as follows: 


1. coSs(pi/6 ) (ans = 0.8660) 
2. tan(45*pi/180 ) (ans = 1.0000) 
3. Sin(pi)+cos(45*pi/180 ) (ans = 0.7071) 


The 
format 


Function 


The format function is used to control how the numeric values are 
displayed in the Command Window. The short format is set by default 


and the numerical results are displayed with 4 digits after the decimal point 
(see the examples above). The Long format produces 15 digits after the 
decimal point. 


Example: 
Calculate 0 = tan(4) and display results in short and Long formats. 
The short format is set by default: 


>> theta=tan(pi/3) 
theta = 
732i 
>> 
And the Long format is tumed on by typing Format long: 
>> theta=tan(pi/3) 
theta = 
W. 732i 


>> format long 
>> theta 


theta = 


1.732050807568877 


Variables 


In MATLAB, a named value is called a variable. MATLAB comes with 
several predefined variables. For example, the name pi refers to the 
mathematical quantity m, which is approximately pl ans = 3.1416 


Note:MATLAB is case-sensitive, which means it distinguishes between 
upper- and lowercase letters (e.g. data, DATA and DaTa are three different 
variables). Command and function names are also case-sensitive. Please 
note that when you use the command-line help, function names are given 
in upper-case letters (e.g., CLEAR) only to emphasize them. Do not use 
upper-case letters when running functions and commands. 


Declaring Variables 


Variables in MATLAB are generally represented as matrix quantities. 
Scalars and vectors are special cases of matrices having size 1x1 (scalar), 
1xn (row vector) or nx1 (column vector). 


Declaration of a Scalar 


The term scalar as used in linear algebra refers to a real number. 
Assignment of scalars in MATLAB is easy, type in the variable name 
followed by = symbol and a number: 


Example: 
a=1 


Comma indo = 


Assignment of a scalar quantity. 


Declaration of a Row Vector 


Elements of a row vector are separated with blanks or commas. 


Example: 
Let's type the following at the command prompt: 
De=" [aise 2 3 4 5] 


'Command Window 


>> b=(1 2 3 4 5] 


Assignment of a row vector quantity. 


We can also use the New Variable button to assign a row vector. In the tool 
strip, select Home > New Variable. This action will create a variable called 
unnamed which is displayed in the workspace. By clicking on the title 
unnamed, we can rename it to something more descriptive. By double- 


clicking on the variable, we can open the Variable Editor and type in the 
values into spreadsheet looking table. 


L Temporaryltems 
» 2008 
» builds 
} Custom Office Templates 


-OmMMandG MIstory 
%-- 03/06/2015 09:... 


v 


Select a file to view details 


Using the New Variable button in the tool strip. 


lolx] 


@ \) $RECYCLE.BIN 
) Temporaryltems 


woont DH wn 


%-- 03/06/2015 09:... 


=y 
So 


clc 


pe —— = 7 | b 
me Ol nd Windo\ 


Select a file to view details 


Assignment of a row vector by using the Variable Editor. 


Declaration of a Column Vector 


Elements of a column vector is ended by a semicolon: 


Example: 
cf) 27374,5,1 


>> c=[1727374757] 


i 
2 
3 
4 
5 


Assignment of a column vector quantity. 


Or by transposing a row vector with the ' operator: 
G =4[1- 273945] 


>> o=[1 2 34 5]* 


c= 


Ob © he be 


Assignment of a column vector quantity by transposing a row vector 
with the ' operator. 


Or by using the Variable Editor: 


Eaigrpa 


@ J) $RECYCLE.BIN 
L Temporaryltems = [1,2,3,4,5] 
® )) 2008 : SE 11;23;4:51 
» builds - 

} Custom Office Templates 


OMMANG MIStory 
%-- 03/06/2015 09:... 
clc 


Details 


Select a file to view details 


Assignment of a column vector quantity by using the Variable Editor. 


Declaration of a Matrix 


Matrices are typed in rows first and separated by semicolons to create 
columns. Consider the examples below: 


Example: 
Let us type in a 2x5 matrix: 
dis [27 46-8 1043757, 9] 


>> d=[2 46610; 135 7 9] 


Assignment of a 2x5 matrix. 


a | ei a 


© ) $RECYCLE.BIN ; 

) Temporaryltems : ' [1,2,3,4,5] 

& § 2008 (1;2;3;4;5] 

@ |) builds § ze [2,4,6,8,10;13,5,7... 
} Custom Office Templates 

Desktop 

& |) Documents 


Select a file to view details 


Assignment of a matrix by using the Variable Editor. 


Example: 

This example is a 5x2 matrix: 

Command Window oo oe 
>> e = [2 43 6 B83 10 123 14 167 18 20] 


2 4 
6 8 
10 12 
14 16 
18 20 


Assignment of a 5x2 matrix. 


Linear Equations 


Systems of linear equations are very important in engineering studies. In the 
course of solving a problem, we often reduce the problem to simultaneous 
equations from which the results are obtained. As you learned earlier, 
MATLAB stands for Matrix Laboratory and has features to handle matrices. 
Using the coefficients of simultaneous linear equations, a matrix can be 
formed to solve a set of simultaneous equations. 


Example: 
Let's solve the following simultaneous equations: 
Equation: 

oy — 1 
Equation: 


22 — oy = 9 


First, we will create a matrix for the left-hand side of the equation using 
the coefficients, namely 1 and 1 for the first and 2 and -5 for the second. 
The matrix looks like this: 

it 

lhe =) 


Equation: 
The above matrix can be entered in the command window by typing A=[ 1 
el |. 
Second, we create a column vector to represent the right-hand side of the 
equation as follows: 

1 

9 


Equation: 

The above column vector can be entered in the command window by 
typing B= [1;9]. 

To solve the simultaneous equation, we will use the left division operator 
and issue the following command: C=A\B. These three steps are illustrated 
below: 


>> A=[1 1; 2 -5] 


—— 

il 1 

2 -5 
>> B= [1;9] 
B= 

1 


>> 


The result C indicating 2 and 1 are the values for x and y, respectively. 


Polynomials 


In the preceding section, we briefly learned about how to use MATLAB to 
solve linear equations. Equally important in engineering problem solving is 
the application of polynomials. Polynomials are functions that are built by 
simply adding together (or subtracting) some power functions. (see 
Wikipedia). 

Equation: 


az’ + br +c=0 
Equation: 
f(x) = az? + br +c 


The coeffcients of a polynominal are entered as a row vector beginning with 
the highest power and including the ones that are equal to 0. 


Example: 
Create a row vector for the following function: 
y = 224+ 32° 4+ 5a? +2410 


Notice that in this example we have 5 terms in the function and therefore 
the row vector will contain 5 elements. p=[2 3 5 1 10] 


Example: 

Create a row vector for the following function: y = 3a* + 4x? — 5 

In this example, coefficients for the terms involving power of 3 and 1 are 
0. The row vector still contains 5 elements as in the previous example but 
this time we will enter two zeros for the coefficients with power of 3 and 1: 
p=[3 0 4 0 -5]. 


The 
polyval 
Function 


We can evaluate a polynomial p for a given value of X using the syntax 
polyval(p, x) where p contains the coefficients of polynomial and x is 
the given number. 


Example: 
Evaluate f(x) at 5. 
Equation: 


f(x) = a+ 2x4 1 


The row vector representing f(x) above is p=[3 2 1]. To evaluate f(x) at 
5, we type in: polyval(p, 5). The following shows the Command 
Window output: 


>> p=[3 2 1] 


>> polyval(p,5) 
ans = 
86 


>> 


The 
roots 
Function 


Consider the following equation: 
Equation: 


ax? + br +c=0 


Probably you have solved this type of equations numerous times. In 
MATLAB, we can use the roots function to find the roots very easily. 


Example: 
Find the roots for the following: 
Equation: 


0.627 + 0.32 — 0.9 =0 


To find the roots, first we enter the coefficients of polynomial in to a row 
vector p with p=[0.6 ©.3 -0©.9] and issue the r=roots(p) 
command. The following shows the command window output: 


>> p=[0.6 0.3 -0.9] 
p = 
0.6000 0.3000 -0.9000 


>> r=roots(p) 


>> 


Splitting a Statement 


You will soon find out that typing long statements in the Command 
Window or in the the Text Editor makes it very hard to read and maintain 
your code. To split a long statement over multiple lines simply enter three 
periods "..." at the end of the line and carry on with your statement on the 
next line. 


Example: 
The following command window output illustrates the use of three periods: 


>> Sin(pi)+cos(45*pi/180) - 
Sin(pi/2)+cos(45*pi/180)+tan(pi/3) 


ans = 
2.1463 


>> Sin(pi)t+cos(45*pi/180)-sin(pi/2)... 
+cos(45*pi/180)+tan(pi/3) 


ans = 
2.1463 


>> 


Comments 


Comments are used to make scripts more "readable". The percent symbol % 
separates the comments from the code. Examine the following examples: 


Example: 

The long statements are split to make it easier to read. However, despite 
the use of descriptive variable names, it is hard to understand what this 
script does, see the following Command Window output: 


t_water=80; 
t_outside=15; 
inner_dia=0.05; 
thickness=0.006; 
Lambda_steel=48; 


AlfaInside=2800; 
AlfaOutside=17; 
thickness_insulation=0.012; 
Lambda_insulation=0.03; 


r_i=inner_dia/2 
r_o=r_it+thickness 
r_i_insulation=r_o 
r_o_insulation=r_i_insulation+thickness_insulation 
AreaInside=2*pi*r_i 
AreaOutside=2*pi*r_o 
AreaOutside_insulated=2*pi*r_o_insulation 
AreaM_pipe=(2*pi*(r_o-r_1))/log(r_o/r_i) 
AreaM_insulation=(2*pi*(r_o_insulation- 
r_i_ insulation) ) 
/log(r_o_insulation/r_i_ insulation) 
TotalResistance=(1/(AlfaInside*AreaInside) )+ 
(thickness/(Lambda_steel*AreaM_pipe) )+ 
(1/7(ALfaOutside*AreaOutside) ) 
TotalResistance_insulated= 
(1/7(AlfaInside*AreaInside) )+ ... 
(thickness/(Lambda_steel*AreaM_pipe) )+ 
(thickness_insulation 
/(Lambda_insulation*AreaM_insulation) )+ 
(1/(AlLfaOutside*AreaOutside_insulated) ) 
Q dot=(t_water-t_outside)/(TotalResistance*1000) 
Q dot_insulated=(t_water - 
t_outside)/(TotalResistance_insulated*1000 ) 
PercentageReducttion=((Q_dot- 
Q dot_insulated)/Q_dot)*100 


Example: 
The following is an edited version of the above including numerous 
comments: 


% Problem 16.06 

% Problem Statement 

% Calculate the percentage reduction in heat loss 
when a layer of hair felt 


% 1S wrapped around the outside surface (see 
problem 16.05) 


Format short 


% Input Values 


t_water=80; % Water temperature [C] 
t_outside=15; % Atmospheric temperature [C] 
inner_dia=0.05; % Inner diameter [m] 
thickness=0.006; % [m] 

Lambda_steel=48; % Thermal conductivity of 
steel [W/mK] 

AlfaInside=2800; % Heat transfer coefficient of 
inside [W/m2K] 

AlfaOutside=17; % Heat transfer coefficient of 


outside [W/m2K] 

% Neglect radiation 

% Additional layer 
thickness_insulation=0.012; % [m] 
Lambda_insulation=0.03; % Thermal 
conductivity of insulation [W/mK] 


% Output Values 

% Q dot=(t_water-t_outside)/TotalResistance 

% TotalResistance=(1/(AlLfaInside*AreaInside) )+ 
(thickness/(Lambda_steel*AreaM) )+ 
(1/(ALfaOutside*AreaOutside) 

% Calculating the unknown terms 
r_i=inner_dia/2 

% Inner radius of pipe [m] 

r_o=r_it+thickness 

% Outer radius of pipe [m] 

r_1_insulation=r_o 

% Inner radius of insulation [m] 
r_o_insulation=r_i_insulation+thickness_insulation 
% Outer radius of pipe [m] 


AreaInside=2*pi*r_i 

AreaOutside=2*pi*r_o 
AreaOutside_insulated=2*pi*r_o_insulation 
AreaM_pipe=(2*pi*(r_o-r_1))/log(r_o/r_i) 
% Logarithmic mean area for pipe 
AreaM_insulation=(2*pi*(r_o_insulation- 
r_i_ insulation) ) 

/log(r_o_insulation/r_i_insulation) % 
Logarithmic mean area for insulation 
TotalResistance=(1/(AlfaInside*AreaInside) )+ 
(thickness/ 

(Lambda_steel*AreaM_pipe) )+ 
(1/7(AlLfaOutside*AreaOutside) ) 
TotalResistance_insulated= 
(1/(AlfaInside*AreaInside) )+(thickness/ 

(Lambda_steel*AreaM_pipe) )+ 
(thickness_insulation/(Lambda_insulation*AreaM_ins 
late LON) ) scan: 

+(1/(AlfaOutside*AreaOutside_insulated) ) 

Q dot=(t_water-t_outside)/(TotalResistance*1000) % 
converting into kW 

Q_ dot_insulated=(t_water - 
t_outside)/(TotalResistance_insulated*1000) % 
converting into kW 

PercentageReducttion=((Q_dot- 

Q dot_insulated)/Q_dot)*100 


Basic Operations 


Command 
sum 

prod 

sqrt 

log10 

log 

max 

min 

mean 


std 


Basic operations. 


Meaning 

Sum of array elements 

Product of array elements 
Square root 

Common logarithm (base 10) 
Natural logarithm 

Maximum elements of array 
Minimum elements of array 
Average or mean value of arrays 


Standard deviation 


Special Characters 


Character 


Meaning 
Assignment 
Prioritize operations 


Construct array 


Character Meaning 


% 


Specify range of array elements 

Row element separator in an array 
Column element separator in an array 
Continue statement to next line 

Decimal point, or structure field separator 


Insert comment line into code 


Special Characters 


Summary of Key Points 


1. 


MATLAB has the common functions found on a scientific calculator 
and can be operated in a similar way, 


. MATLAB can store values in variables. Variables are case sensitive 


and some variables are reserved by MATLAB (e.g. Pi stores 3.1416), 


. Variable Editor can be used to enter or manipulate matrices, 
. The coefficients of simultaneous linear equations and polynomials are 


used to form a row vector. MATLAB then can be used to solve the 
equations, 


. The format function is used to control the number of digits 
displayed, 
. Three periods "..." at the end of the line is used to split a long 


statement over multiple lines, 


. The percent symbol % separates the comments from the code, 


anything following % symbol is ignored by MATLAB. 


Problem Set 
Problem Set for MATLAB Essentials 


Determine the value of each of the following. 
Exercise: 


Problem: 6 x 7 + 42 — 24 
Solution: 
>> (6*7)+442-244 (ans = 42) 


Exercise: 


, oe g4°5_52 
Problem: Bost + Bast 
Solution: 


>> ((342+243)/(445-544) )+ 
((sqrt (64) -542)/(445+546+748) ) (ans = 0.0426) 


Exercise: 


Problem: log 107 + 10° 


Solution: 
>> 1log10(1042)+105 (ans = 100002) 
Exercise: 


Problem: e? + 2° — In (e?) 


Solution: 


>> exp(2)+243-log(exp(2) ) (ans = 13.3891) 


Exercise: 


Problem: sin(27) + cos (4) 


Solution: 
>> Sin(2*pi)+cos(pi/4) (ans = 0.7071) 


Exercise: 


Problem: tan( 4) + cos(270) + sin(270) + cos(4) 


Solution: 


>> 
tan(pi73 )+co0s(2/0* p1/180 )+sin( 2/70" pi7/180)+cos(p 
1/3) (ans = 1.2321) 


Exercise: 


Solve the following system of equations: 
22 -+4y=1 
Problem: z + 5y = 2 


Solution: 


SS A=|2 Ao Ss as aa SS Bat aa 
>> Solution=A\B Solution = -0.5000 0.5000 


Exercise: 


Evaluate y at 5. 
Problem: y = 4x4 + 3x7 — x 


Solution: 


>> p=[4 0 3 -1 0] p=40 3 -1 0 >> 
polyval(p,5) ans = 2570 >> 


Exercise: 
Problem: 
Given below is Load-Gage Length data for a type 304 stainless steel 


that underwent a tensile test. Original specimen diameter is 12.7 mm. 
[footnote | 


Load [N] Gage Length [mm] 
0.000 50.8000 
4890 50.8102 
9779 50.8203 
14670 50.8305 
19560 50.8406 
24450 50.8508 
27620 50.8610 
29390 50.8711 
32680 50.9016 
33950 50.9270 
34580 50.9524 


35220 50.9778 


Load [N] Gage Length [mm] 


35720 51.0032 

40540 51.816 

48390 53.340 

59030 55.880 

65870 58.420 

69420 60.960 

69670 (maximum) 61.468 

68150 63.500 

60810 (fracture) 66.040 (after fracture) 

= £, where P is the load [N] on the sample with an original cross- 


sectional area A [m7] and the engineering strain is defined as € = oar 


where Al is the change in length and J is the initial length. 


Compute the stress and strain values for each of the measurements 
obtained in the tensile test. Data available for download. 
Introduction to Materials Science for Engineers by J. F. Shackelford, 
Macmillan Publishing Company. © 1985, (p.304) 


Solution: 


First, we need to enter the data sets. Because it is rather a large table, 
using Variable Editor is more convenient. See the figures below: 


7] | ED No 


~“$-— 30/10/2011 12:03 --% 


Ee A 


<21x1 double> 
<21x1 double> 


~~$-— 30/10/2011 12:03 --% 


Extension length in mm. 


Next, we will calculate the cross-sectional area. Area=pi/4* 
(©.012742) Area = 1.2668e-004 


Now, we can find the Stress values with the following, note that we are 
obtaining results in MPa: Sigma=(Load_N./Area) *104(-6) 


Sigma = 0 38.6022 77.1964 115.8065 154.4086 
1930108 216.0351) 232.0076 257.9792 263.0047 
Ziz2.9700) 270.0502 201,97 735 320), 0269) 361, 9955 
465.9888 519.9844 548.0085 549.9820 537.9830 
480.0403 


For strain calculation, we will first find the change in length: 
Delta_L=Length_mm-50.800 Delta_L = 0 0.0102 
0.0208 0.0305 ©:0406 ©.0508 0.0610 0.0711 
O2016 O72 7/0 021524 0.11/73 022032 1.0160 
2.5400 5.0800 7.6200 10.1600 10.6680 12.7000 
15.2400 


Now we can determine Strain with the following: 
Epsilon=Delta_L./50.800 Epsilon = 0 0.0002 
0.0004 0.0006 0.0008 0.0010 0.0012 0.0014 
6.0020 ©, 0025 0.0030 0.0035 0.0040 0.0200 
020506 07; 1000) €21500 ©. 2000 0.2100 052500 
0.3000 


The final results can be tabulated as foolows: [Sigma Epsilon] 


ans = 0 0 38.6022 0.0002 77.1964 0.0004 

115.8065 0.0006 154.4086 0.0008 193.0108 0.0010 
213, 0351.0 ,0012 (232, 0076 .0 0014 1257.97 92 60,0020 
260,0047 ©0025 272,970 0, 0030 278.0302 ©,0035 
201,97 73 (0).0040 32070269 0.0200) 381.9955 0.0500 
465.9888 0.1000 519.9844 0.1500 548.0085 0.2000 
949.9820 0.2100 537.9830 0.2500 480.0403 0.3000 


Plotting in MATLAB 
Plotting basics. 


2-0 Plo+s 
2D Plots 
Annotation of Plot. 


& 


www. hetemeel.com 


A picture is worth a thousand words, particularly visual representation of 

data in engineering is very useful. MATLAB has powerful graphics tools 

and there is a very helpful section devoted to graphics in MATLAB Help: 
Graphics. Students are encouraged to study that section; what follows is a 
brief summary of the main plotting features. 


Two-Dimensional Plots 


The 
plot 
Statement 


Probably the most common method for creating a plot is by issuing 
plot(x, y) statement where function y is plotted against x. 


Example: 
Type in the following statement at the MATLAB prompt: 


X=[-pi:.1:pi]; y=sin(x); plot(x,y),; 
After we executed the statement above, a plot named Figure1 is generated: 


| 
0.8 
0.6 


0.4 


Graph of sin(x) 


Having variables assigned in the Workspace, x and y=sin(x) in our case, we 
can also select x and y, and right click on the selected variables. This opens 
a menu from which we choose plot(x,y). See the figure below. 


<1x63 double> 


i ~ = 
nenamMe 


EGit Value 


= 


Creating a plot from 
Workspace. 


Annotating Plots 


Graphs without labels are incomplete and labeling elements such as plot 
title, labels for x and y axes, and legend should be included. Using up 


arrow, recall the statement above and add the annotation commands as 
shown below. 


X=[-pi:.1:p1];y=sin(x);plot(x,y);title('Graph of 
y=sin(x)');xlabel('x');ylabel('sin(x)')jgrid on 


Run the file and compare your result with the first one. 


Graph of y=sin(x) 


sin(x) 


Graph of sin(x) with Labels. 


Note: Type in the following at the MATLAB prompt and learn additional 
commands to annotate plots: 


help gtext 
help legend 
help zlabel 


Superimposed Plots 


If you want to merge data from two graphs, rather than create a new graph 
from scratch, you can superimpose the two using a simple trick: 


% This script generates sin(x) and cos(x) plot on 
the same graph 
% initialize variables 


X=[-p1l:.1:pli]; %create a row vector from -pi 
to +pi with .1 increments 

yO=sin(x); %cCalculate sine value for each 
x 

yi=cos(x); %Calculate cosine value for 
each x 


% Plot sin(x) and cos(x) on the same graph 
plot(x,y0,x,y1); 

title('Graph of sin(x) and cos(x)'); %Title of 
graph 

Xlabel('x'); %Label of x axis 
ylabel('sin(x), cos(x)'); %Label of y axis 
legend('sin(x)', 'cos(x)'); %Insert legend in the 
same order as yO and yi calculated 

grid on %Graph grid is turned 


Graph of sin(x) and cos(x) 


sin(x), cost) 


Graph of sin(x) and cos(x) in the same plot with labels and 
legend. 


Multiple Plots in a Figure 


Multiple plots in a single figure can be generated with Subp1ot in the 
Command Window. However, this time we will use the built-in Plot Tools. 
Before we initialize that tool set, let us create the necessary variables using 
the following script: 


% This script generates sin(x) and cos(x) 
variables 

clc %Clears command window 
Clear all %Clears the variable space 


Close all %Closes all figures 
X1=[-2*pi:.1:2*pi]; %Creates a row vector from 
-2*pi to 2*pi with .1 increments 


Y1=sin(X1); %Calculates sine value for 
each x 

Y2=cos(X1); %Calculates cosine value for 
each x 

Y3=Y1+Y2; %Calculates sin(x)+cos(x) 
Y4=Y1-Y2; %Calculates sin(x)-cos(x) 


Note that the above script clears the command window and variable 
workspace. It also closes any open Figures. After running the script, we will 
have X1, Y1, Y2, Y3 and Y4 loaded in the workspace. Next, select File > 
New > Figure, a new Figure window will open. Click "Show Plot Tools and 
Dock Figure" on the tool bar. 


a tan (pi/3 
ale 


Plot Tools 


Under New Subplots > 2D Axes, select four vertical boxes that will create 
four subplots in one figure. Also notice, the five variables we created earlier 
are listed under Variables. 


Creating four sub plots. 


After the subplots have been created, select the first supblot and click on 
"Add Data". In the dialog box, set X Data Source to X1 and Y Data Source 
to Y1. Repeat this step for the remaining subplots paying attention to Y 
Data Source (Y2, Y3 and Y4 need to be selected in the subsequent steps 
while X1 is always the X Data Source). 


Adding data to axes. 


Next, select the first item in "Plot Browser" and activate the "Property 
Editor". Fill out the fields as shown in the figure below. Repeat this step for 
all subplots. 


Mw —YlvsX1 


M Axes (no title) 


Using "Property Editor". 


Save the figure as Sinxcosx. fig in the current directory. 


F Graph of sin(x) 
— Yl vs X1 


Graph of cos(x) 
— Y2 vs X1 


Graph of sin(x)+cos(x) 
M —Y3 vs X1 


M Graph of sin(x)-cos(x) 
MW —Y4vsX1 


The four subplots generated with "Plot Tools". 


The four subplots in a single figure. 


Three-Dimensional Plots 


3D plots can be generated from the Command Window as well as by GUI 
alternatives. This time, we will go back to the Command Window. 


The 
plot3 
Statement 


With the X1,Y1,Y2 and Y2 variables still in the workspace, type in 
plot3(X1, Y1, Y2) at the MATLAB prompt. A figure will be generated, 
click "Show Plot Tools and Dock Figure". 


A raw 3D figure is generated with plots. 


Use the property editor to make the following changes. 


3D Property Editor. 


The final result should look like this: 


Figures - Figure 1 


NGM e/k|(AAT9Be : 


3D graph of x, sin(x), cos(x) 


Use help or doc commands to learn more about 3D plots, for example, 
image(x), surf(x) andmesh(x). 


Quiver or Velocity Plots 


To plot vectors, it is useful to draw arrows so that the direction of the arrow 
points the direction of the vector and the length of the arrow is vector’s 
magnitude. However the standard plot function is not suitable for this 
purpose. Fortunately, MATLAB has quiver function appropriately named 
to plot arrows. quiver (xX, Y,U,V) plots vectors as arrows at the 
coordinates (x,y) with components (u,v). The matrices x, y, u, and v must 


all be the same size and contain corresponding position and velocity 
components. 


Example: 


Calculate the magnitude of forces OA, OB and the resultant R of OA and 
OB shown below. Plot all three forces on x-y Cartesian coordinate 
system|[footnote]. 

Ay 


A (600,320) 


0 (0,0) x 


B (-200,-480) 


Quiver Plot. 


Applied Engineering Mechanics by A. Jensen, H. Chenoweth McGraw-Hill 
Ryerson Limited © 1972, (p. 15) 


% Preparation 


clear % removes all variables from the current 
workspace, 


% releasing them from system memory. 
cle % Clears all input and output from the 
Command Window display, 

% giving you a "clean screen." 
% Input and Computation 
OA=[600 320]; % Force 1 
magOA=sqrt(sum(0A.42)); 
OB=[-200 -480]; % Force 2 
magOB=sqrt(sum(0B.42)); 
OC=OA+0B; % The resultant of OA and OB 
magOC=sqrt(sum(0OC.42)); % The magnitude of 
resultant force OC 
angleMag=atan(0C(2)/0C(1))*180/pi; % angle of OC 
in degrees 
% Output 
disp(' ') % Display blank line 
stri= ['The magnitude of the resultant force is ', 
num2str(magOC), ' N.']; 
disp(stri); 
str2= ['The angle of the resultant force is ', 
num2str(angleMag), ' degrees.']; 
disp(str2); 
% Plot Preparation 
starts = zeros(3,2); % Origin for all 3 forces, 
3x2 "zero" matrix 
ends = [0A;0B;0C]; % End point for all 3 forces 
vectors = horzcat(starts,ends); % Concatenate 
arrays horizontally 
% Plot Forces on x-y Cartesian Coordinate System 
% The following MATLAB function plots vectors as 
arrows 
% at the coordinates specified in each 
corresponding 
% pair of elements in x and y. 
quiver( vectors( :,1 ), vectors( :,2 ), vectors( 
po Jy VECTORS( <,4 2))- 
axis equal 


grid 

title('Forces on x-y Cartesian Coordinate System' ) 
Xlabel('x') % x-axis label 

ylabel('y') % y-axis label 

view(2) % setting view to 2-D 


Forces on x-y Cartesian Coordinate System 
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Output of quiver function. 
Example: 


Write an interactive script to calculate the resultant R of forces F1, F2 and 
F3 shown below and plot all four forces on x-y Cartesian coordinate 
system|footnote]. 


F2=200 N 
at 150 deg 


F3=330 N 
at 270 deg 


F1=170N 
at 60 deg 


An example for quiver3 plot. 


v x 


Applied Engineering Mechanics by A. Jensen, H. Chenoweth McGraw-Hill 
Ryerson Limited © 1972, (p. 15) 


clear 
ele 


disp('This script computes the resultant of three 
forces on x-y Cartesian coordinate system. ') 


fi=input('Enter the 
‘); 
thetal=input('Enter 
deg: '); 
f2=input('Enter the 
Ni); 
theta2=input('Enter 
deg: ‘); 
f3=input('Enter the 
‘); 

theta3=input (‘Enter 
deg: ‘); 


magnitude 
the angle 
magnitude 
the angle 
magnitude 


the angle 


of 


of 


of 


of 


of 


of 


first force in N: 


first force in 


second force in 


second force in 


third force in N: 


third force in 


x1=f1*cos(thetai*pi/180); % The components of 
force 

yi=fi*sin(theta1*pi/180); % The components of 
force 

F1i=[x1 y1]; % Force 1 

x2=f2*cos(theta2*pi/180); % The components of 
force 

y2=f2*sin(theta2*pi/180); % The components of 
force 

F2=[x2 y2]; % Force 2 

X3=fF3*cos(theta3*pi/180); % The components of 
force 

y3=f3*sin(theta3*pi/180); % The components of 
force 

F3=[x3 y3]; % Force 3 

R=F1+F2+F3; % The resultant of Fi, F2 and F3 
magR=sqrt(sum(R.42)); % The magnitude of resultant 
force R 

angle=atan(R(2)/R(1))*180/pi; % Angle of R in 
degrees 

disp(' ') % Display blank line 
stri= ['The magnitude of the resultant force is ', 
num2str(magR), ' N.']; 

disp(stri); 

str2= ['The angle of the resultant force is ', 
num2str(angle), ' degrees.']; 

disp(str2); 

starts = zeros(4,3); 

ends = [F1;F2;F3;R]; 

ends(3,3)=0; % inputs Os for zZ components, making 
it 3D 

vectors = horzcat(starts,ends); % Concatenate 
arrays horizontally 

quiver3( vectors( :,1 ), vectors( :,2 ), vectors( 
:,3 ), vectors( :,4 ), vectors( :,5 ), vectors( 
:,6 )); % A three-dimensional quiver plot displays 
vectors with components (u,v,w) at the points 


(x,y,Z), where u, v, w, X, y, and z all have real 
(non-complex) values. 

axis equal 

title('Forces on x-y Cartesian coordinate system’ ) 
Xlabel('x') % x-axis label 

ylabel('y') % y-axis label 

view(2) 


Forces on x-y Cartesian coordinate system 
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Output of quiver 3 function. 


Summary of Key Points 


1.plot(x, y) and plot3(X1, Y1, Y2) statements create 2- and 3- 
D graphs respectively, 

2. Plots at minimum should contain the following elements: title, 
xlabel, ylabel and legend, 

3. Annotated plots can be easily generated with GUI Plot Tools, 

4. quiver and quivers plots are useful for making vector diagrams. 


Problem Set 
Problem Set for Graphing with MATLAB 
Exercise: 


Problem: 


Plot y = a + bz, using the specified coefficients and ranges (use 


increments of 0.1): 


aa=2,b=0,3ior0 <2 5 
b.a= 3,b= 0 for0 < 2 < 10 
6a=4,b=—0.3 ford <2 < 15 


Solution: 


a.a=2; b=.3; x=[0:.1:5]; y=atb*x 
DIGE(X, ¥)> eLele( Graph an 
y=atbx'),xlabel('x'),ylabel('y'),grid 

b.a=3; b=.0; x=[0:.1:10]; y=at+b*x 
DEGE(G;Y )), ELele(? Graph Om 
y=atbx'),xlabel('x'),ylabel('y'),grid 

c.a=2; b=.3; x=[0:.1:5]; y=atb*x 
plot(x,y),title('’Graph of 
y=atbx'),xlabel('x'),ylabel('y'),grid 


Exercise: 


Problem: 


Plot the following functions, using increments of 0.01 and a = 6, 


b=0.8,0 =< 2 < 5: 


ay=at+2° 
b.y= az? 
c. y = asin(z) 


Solution: 
a. a=6; b=.8; x=[0:.01:5]; y=atx.Ab; 
DIOE(X, y), titlel Graph of 
y=atxAb'),xlabel('x'), ylabel('y'),grid 


Graph of y=atx? 


b.a=6; b=.8; x=[0:.01:5]; y=a*x.Ab; 
plot(x,y),title('Graph of 
y=axAb'),xlabel('x'),ylabel('y'),grid 


Graph of y=ax” 


25 


c.a=6; x=[0:.01:5]; y=a*sin(x); 
plot(x,y),title('Graph of 
y=a*sin(x)'),xlabel('x'),ylabel('y'),grid 


Graph of y=a*sin(x) 
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Exercise: 
Problem: 
Plot function y = sin) for =i; <x < 107 using increments of =75 


Solution: 


= p10 p77 100 Opi Vy = sini x 
plot(x,y),title('Graph of 
y=sin(x)/x'),xlabel('x'),ylabel('y'),grid 


sin(z) 


Graph of y = 


Exercise: 


Problem: 


Data collected from Boyle's Law experiment are as follows: (Data 
available for download.) 


Volume [cmA3] Pressure [Pa] 


7.34 100330 
7.24 102200 
7.14 103930 
7.04 105270 
6.89 107400 
6.84 108470 
6.79 109400 
6.69 111140 
6.64 112200 


Plot a graph of Pressure vs Volume, annotate your graph. 
Solution: 


Pressure= 

(186330, 102260) 103930, 105270, 1107400), 108476, 1094 
00,111140,112200]; Volume= 

(Wone4nt 24,7. 14-7 04, 6.696.084, 6.79) 616976.04)- 
plot(Volume, Pressure), title('Pressure Volume 
Graph'),xlabel('Volume'),ylabel( 'Pressure'),gri 
d 


«10° Pressure Yolume Graph 


Pressure 


Volume 


Exercise: 


Problem: 


The original data collected from Boyle's [footnote] experiment are as 
follows: (Data available for download.) 


Volume [tube-inches] Pressure [inches-Hg] 
12 29.125 


10 35.000 


Volume [tube-inches] Pressure [inches-Hg] 


8 43.688 
6 98.250 
rs) 70.000 
4 87.375 
3 116.500 


Plot a graph of Pressure vs Volume, annotate your graph. 
Introduction to Engineering: Modeling and Problem Solving by J. B. 
Brockman, John Wiley and Sons, Inc. © 2009, (p.246) 


Solution: 


ee P=)29 7125, 357432680, 58.25, 70, 872275, l1b25] 
pe V=/(12,10,8,6,5,4,317 2 > 

plot(vV,P),title( ‘Pressure Volume 
Graph'),xlabel( 'Volume'),ylabel( 'Pressure'),gri 
d 


Pressure Volume Graph 
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Exercise: 


Problem: Display the two plots created earlier in one plot. 


Solution: 


Original Experiment 


Pressure 


Volume 


x 10 Experimental Data 


Pressure 


Exercise: 


Problem: 


A tensile test of SAE 1020 steel produced the data below (Data 
available for download.) [footnote] experiment are as follows: 


Extension [mm] Load [kN] 


0.00 0.0 


Extension [mm] Load [KN] 


0.09 1.9 
0.31 6.1 
0.47 9.4 
2.13 11.0 
9.05 11.7 
10.50 12.0 
16.50 11.9 
23.70 10.7 
27.70 9.3 
34.50 8.1 


Plot a graph of Load vs Extension, annotate your graph. 
Introduction to Materials Science for Engineers | Instructor's Manual 
by J. F. Shackelford, Macmillan Publishing Company. © 1992, (p.440) 


Solution: 


Extension= 

(Oe 007070070) 31) O04, 271375,05, 10 50, 1G, 50, 2347 
0, 27 270734250; Eoad= 

ROL Opt or6r 1) S47 1st Oey 1) G7, one 3 
-1]; plot(Extension, Load),title('Load versus 
Extension 

Curve'),xlabel( 'Extension'), ylabel('Load'), grid 


Load versus Extension Curve 
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Exercise: 


Problem: 


Given below is Stress-Strain data for a type 304 stainless steel. 
[footnote] experiment are as follows: (Data available for download.) 


Stress [MPa] Strain [mm/mm] 
0.0 0.0000 
38.6 


0.0002 


Stress [MPa] Strain [mm/mm|] 


T72 0.0004 
115.8 0.0006 
154.4 0.0008 
193.0 0.0010 
218.0 0.0012 
232.0 0.0014 
258.0 0.0020 
268.0 0.0025 
273.0 0.0030 
278.0 0.0035 
282.0 0.0040 
320.0 0.0200 
382.0 0.0500 
466.0 0.1000 
520.0 0.1500 
548.0 0.2000 


5930.0 0.2100 


Stress [MPa] Strain [mm/mm] 
538.0 0.2500 


480.0 0.3000 


Plot a graph of Stress vs Strain, annotate your graph. 
Introduction to Materials Science for Engineers by J. F. Shackelford, 
Macmillan Publishing Company. © 1985, (p.304) 


Solution: 


The data can be entered using Variable Editor: 


(4 Variable Editor - Stress ra oe Workspace 
aa | & 2 | S| é ~ | te | [ees =] | FO no vaiid pl.» io | ax 


EH Stress <21x1 double> 


115.8065] 
154.4086] 
193.0108] 
218.0351] 
232.0076] 
257.9792| 


Then execute the following: 
plot(Strain, Stress), title('Stress versus Strain 
Curve'),xlabel('Strain [mm/mm]'), ylabel('Stress 
[mPa]'),grid 


Stress [mPa] 
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Stress versus Strain Curve 
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Writing Scripts to Solve Problems 
Basic MATLAB programming concepts are presented to demonstrate how 
to create, save and execute script files. 


Scr ipt files 


Interac Tivit y 
Loops 


Sty le Guidelines 


www, hetemeel.com 


MATLAB provides scripting and automation tools that can simplify 
repetitive computational tasks. For example, a series of commands executed 
in a MATLAB session to solve a problem can be saved in a script file called 
an m-file. An m-file can be executed from the command line by typing the 
name of the file or by pressing the run button in the built-in text editor tool 
bar. 


Script Files 


A script is a file containing a sequence of MATLAB statements. Script files 
have a filename extension of .m. By typing the filename at the command 
prompt, we can run the script and obtain results in the command window. 


ide H:\MATLAB\M-Files_HeatTransfer 
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ThermalConductivity.m (MATLABS 


Stainless Steel 


Number of m-files are displayed in the Current Folder sub-window. 


A sample m-file named ThermalConductivity .m is displayed in Text 
Editor below. Note the triangle (in green) run button in the tool bar, pressing 
this button executes the script in the command window. 


ATLAB\M-Files_HeatTransfer\ ThermalConductivity.m 


File Edit Text Go Cell Tools Debug Desktop Window Help 3 || x 


PO Gie|se@arelou-|\mevh &)- 22 B® BT BB | stack| ease ~]| & BoOeaes|o 
7 dOpen fie Cui-O)s fry x [a /@ 


% Stainless Steel 

lambda_StainlessSteel=14.4; % Thermal Conductivity of Staninless Steel [W/mK] 
deltaZ_1=1; % Distance between the first set of thermocouples [cm] 
deltaZ_2=1; % Distance between the second set of thermocouples [cm] 
deltaZ_SS_average=(deltaZ_l+deltaZ_2)/2; % Average Distance [cm] 

t1i=440; % Temperature 1 [C] 

t2=400; % Temperature 2 [C] 

t3=360; % Temperature 3 [C] 

deltaT_12=ti-t2; % Temperature difference between thermocouples 1 and 2 [C] 
deltaT_23=t2-t3; % Temperature difference between thermocouples 2 and 3 [C] 
deltaT_SS_average=(deltaT_12+deltaT_23)/2; % Average temperature difference for Stainless Steel [C] 


= Aluminum 

deltaZ_3=3; % Distance between the third set of thermocouples [cm] 
deltaZ_4=3; % Distance between the fourth set of thermocouples [cm] 
deltaZ_Al_ average=(deltaZ 3+deltaZ_4)/2; % Average Distance [cm] 

£4=270: % Temperature 4 [C] 

£S=261.47 % Temperature 5 [C] 

t6=252.7;7 % Temperature 6 [C] 

deltaT_45=t4-t5; % Temperature difference between thermocouples 4 and 5 [C] 
deltaT_Sé6=t5-té; % Temperature difference between thermocouples 5 and 6 [C] 
deltaT_Al_average=(deltaT_45+deltaT_56)/2; % Average temperature difference for Aluminum [C] 


% Calculation for Thermal Conductivity of Aluminum [W/mK] 
lambda_Aluminum=lambda_StainlessSteel* (deltaT_SS_average/deltaT_Al_average)*(deltaZ_Al_average/deltaZ_SS_average) 


script In 1 col [ove , 


The content of ThermalConductivity .m file is displayed in Text 
Editor. 


Now let us see how an in-file is created and executed. 


Example: 

A cylindrical acetylene bottle with a radius r=0.3 m has a hemispherical 
top. The height of the cylindrical part is h=1.5 m. Write a simple script to 
calculate the volume of the acetylene bottle. 

To solve this problem, we will first apply the volume of cylinder equation. 
Using the volume of sphere equation, we will calculate the volume of 
hemisphere. The total volume of the acetylene bottle is found with the sum 
of volumes equation. 

Equation: 


2 
Veylinder = ar it 


Equation: 
a 3 
Vsphere = uh 
3 
Equation: 
23 
Viop = ae 
Equation: 


Vacetylene bottle = Vejinder ug Viop 


To write the script, we will use the built-in text editor. From the menu bar 
select File > New > Script. The text editor window will open in a separate 
window. First save this fileas AcetyleneBottle.m. In that window 
type the following code paying attention to the use of percentage and 
semicolon symbols to comment out the lines and suppress the output, 
respectively. 


% This script computes the volume of an acetylene 
bottle with a radius r=0.3 m, 

% a hemispherical top and a height of cylindrical 
part h=1.5 m. 


r=0.3; % Radius [m] 
h=1.5; % Height [m] 
Vol_top=(2*p1i*r43)/3; % Calculating the 
volume of hemispherical top [m3] 
Vol_cyl=pi*r42*h; % Calculating the 


volume of cylindrical bottom [m3] 
Vol_total=Vol_top+Vol_cyl % Calculating the 
total volume of acetylene bottle [m3] 


% This script computes the volume of an acetylene bottle with a radius r=0.3 m, a hemispherical top 
% and a height of cylindrical part h=1.5 m. 


r=0.3; % Radius [m] 


h=1.5; % Height [m] 
Vol_top=(2*pi*r*3) /3; 
Vol_cyl=pi*r“2*h; 

Vol_tot a1=Vo 1_top+Vol_cyl 


% Calculating the volume of hemispherical top [m3] 
% Calculating the volume of cylinderical bottom [m3] 
% Calculating the total volume of acetylene bottle [m3] 


Script created with the built-in text editor. 


After running the script by pressing the green button in the Text Editor tool 
bar, the output is displayed in the command window as shown below. 
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AcetyleneBottle.m (MATLAB Script) ¥ 


This script computes the 
volume of an acetylene bottle 
with a radius r=0.3 m, a 
hemispherical top 


4 start 


The MATLAB output in the command window. 


The 

input 

Function 

Notice that the script we have created above is not interactive and computes 
the total volume only for the variables defined in the m-file. To make this 
script interactive we will make some changes to the existing 
AcetyleneBottle.m by adding input function and save it as 
AcetyleneBottleInteractive.m. 


The syntax for input is as follows: 


userResponse = input('prompt' ) 


Example: 

Now, let's incorporate the Lnput command in 
AcetyleneBottleInteractive.m as shown below and the 
subsequent figure: 


% This script computes the volume of an acetylene 
bottle 
% user is prompted to enter 

% a radius r for a hemispherical top 

% a height h for a cylindrical part 
r=input('Enter the radius of acetylene bottle in 
meters '); 
h=input('Enter the height of cylindrical part of 
acetylene bottle in meters '); 


Vol_top=(2*p1i*r43)/3; % Calculating the 
volume of hemispherical top [m3] 
Vol_cyl=pi*r42*h; % Calculating the 
volume of cylindrical bottom [m3] 
Vol_total=Vol_top+Vol_cyl % Calculating the 


total volume of acetylene bottle [m3] 


% This script computes the volume of an acetylene bottle 
% user is prompted to enter 

% a radius r for a hemispherical top 

% a height h for a cylindrical part 
r=input ('Enter the radius of acetylene bottle in meters '); 
h=input('Enter the height of cylinderical part of acetylene bottle in meters '); 
Vol_top=(2*pi*r*3) /3; % Calculating the volume of hemispherical top [m3] 
Vol_cyl=pi*r*2*h; % Calculating the volume of cylinderical bottom [m3] 
Vol_total=Vol_top+Vol_cyl % Calculating the total volume of acetylene bottle [m3] 

| 


Interactive script that computes the volume of acetylene cylinder. 


The command window upon run will be as follows, note that user keys in 
the radius and height values and the same input values result in the same 
numerical answer as in example which proves that the computation is 
correct. 
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ad ree 
AcetyleneBottle.m (MATLAB Script) ae 
This script computes the he ai 
volume of an acetylene bottle 0.4807 
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4 Star| 


The same numerical result is obtained through interactive script. 


The 
disp 


Function 


As you might have noticed, the output of our script is not displayed in a 
well-formatted fashion. Using disp, we can control how text or arrays are 
displayed in the command window. For example, to display a text string on 
the screen, type in disp('Hello world! ' ). This command will return 
our friendly greeting as follows: Hello world! 


disp(variable) can be used to display only the value of a variable. To 
demonstrate this, issue the following command in the command window: 


b= [12345] 


We have created a row vector with 5 elements. The following is displayed 
in the command window: 


>> b= [12345] 
b= 
1 2 3 4 5 


Now if we type in disp(b) and press enter, the variable name will not be 
displayed but its value will be printed on the screen: 


>> disp(b) 
al 2 | 4 5 


The following example demonstrates the usage of disp function. 


Example: 

Now, let's open ACetyleneBottleInteractive.m file and modify 
it by using the disp command. First save the file as 
AcetyleneBottleInteractiveDisp.m, so that we don't 
accidentally introduce errors to a working file and also we can easily find 
this particular file that utilizes the disp command in the future. The new 
file should contain the code below: 


% This script computes the volume of an acetylene 
bottle 
% user is prompted to enter 
% a radius r for a hemispherical top 


% a height h for a cylindrical part 
cle % Clear screen 
disp('This script computes the volume of an 
acetylene bottle' ) 
r=input('Enter the radius of acetylene bottle in 
meters '); 
h=input('Enter the height of cylindrical part of 
acetylene bottle in meters '); 


Vol_top=(2*p1i*r43)/3; % Calculating the 
volume of hemispherical top [m3] 
Vol_cyl=pi*r42*h; % Calculating the 


volume of cylindrical bottom [m3] 
Vol_total=Vol_top+Vol_cyl; % Calculating the 
total volume of acetylene bottle [m3] 

disp(' ') % Display blank line 
disp('The volume of the acetylene bottle is') % 
Display text 

disp(Vol_total) % Display variable 


Your screen output should look similar to the one below: 


This script computes the volume of an acetylene 
bottle 

Enter the radius of acetylene bottle in meters .3 
Enter the height of cylindrical part of acetylene 
bottle in meters 1.5 


The volume of the acetylene bottle is 
0.4807 


The 


num2str 


Function 


The num2str function allows us to convert a number to a text string. 
Basic syntax is Str = num2str(A) where variable A is converted to a 
text and stored in str. Let's see how it works in 
AcetyleneBottleInteractiveDisp.m. Remember to save the file 
with a different name before editing it, for example, 
AcetyleneBottleInteractiveDisp1.m. 


Example: 
Add the following line of code to your file: 


str = ['The volume of the acetylene bottle is ', 
num2str(Vol_total), ' cubic meters.']; 


Notice that the three arguments in Str are separated with commas. The 
first argument is a simple text that is contained in''. The second argument 
is where the number to string conversion take place. And finally the third 
argument is also a simple text that completes the sentence displayed on the 
screen. Using semicolon at the end of the line suppresses the output. In the 
next line of our script, we will call str with disp(str);. 
AcetyleneBottleInteractiveDisp1.m file should look like this: 


% This script computes the volume of an acetylene 
bottle 
% user is prompted to enter 

% a radius r for a hemispherical top 

% a height h for a cylindrical part 
cle % Clear screen 
disp('This script computes the volume of an 
acetylene bottle:') 
disp(' ') % Display blank line 
r=input('Enter the radius of acetylene bottle in 
meters '); 
h=input('Enter the height of cylindrical part of 


acetylene bottle in meters '); 


Vol_top=(2*p1i*r43)/3; % Calculating the 
volume of hemispherical top [m3] 
Vol_cyl=pi*r42*h; % Calculating the 


volume of cylindrical bottom [m3] 
Vol_total=Vol_top+Vol_cyl; % Calculating the 
total volume of acetylene bottle [m3] 

disp(' ') % Display blank line 
str = ['The volume of the acetylene bottle is ', 
num2str(Vol_total), ' cubic meters.']; 

disp(str); 


Running the script should produce the following: 


This script computes the volume of an acetylene 
bottle: 


Enter the radius of acetylene bottle in meters .3 
Enter the height of cylindrical part of acetylene 
bottle in meters 1.5 


The volume of the acetylene bottle is 0.48066 
Cubic meters. 


The 
fopen 
and 
fclose 


Functions 


The first command is used to open or create a file. The basic syntax for 
fopen is as follows: 


fid = fopen(filename, permission) 


For example, fo = fopen('output.txt', 'w'); opens or creates 
a new file named Output. txt and sets the permission for writing. If the 
file already exists, it discards the existing contents. 


fclose command is used to close a file. For example, if we type in 
fclose(fo);, we close the file that was created above. 


The 
fprintf 
Function 


fprintf function writes formatted data to the computer monitor or a file. 
This command can be used to save the results of a calculation to a file. To 
do this, first we create or open an output file with Fopen, second we issue 
the fprintf command and then we close the output file with fclose. 


The simplified syntax for fprintf is as follows: 


fprintf=(fid, format, variablei, variable 2, ...) 


Example: 
Add the following lines to your .m file: 


fo = fopen('output.txt', 'w'); 

fprintf(fo,'The radius of acetylene bottle: %g 
meters \n', r); 

fprintf(fo, 'The height of cylindrical part of 
acetylene bottle: %g meters \n', h); 


fprintf(fo,'The volume of the acetylene bottle: %g 
Cubic meters. \n', Vol_total); 
fclose(fo); 


Here, we first create the output.txt file that will contain the following 
three variables r, h and Vol_total. In the fo output file, the variables 
are formated with %g which automatically uses the shortest display format. 
You can also use %1 or %d for integers and %e for scientific notation. In 
our script above, the \n (newline) moves the cursor to the next line. 
Naming the new .m file as 
AcetyleneBottleInteractiveOutput .m, it should look like this: 


% This script computes the volume of an acetylene 
bottle 
% user is prompted to enter 

% a radius r for a hemispherical top 

% a height h for a cylindrical part 
(ond res % Clear screen 
disp('This script computes the volume of an 
acetylene bottle:') 
disp(' ') % Display blank line 
r=input('Enter the radius of acetylene bottle in 
meters '); 
h=input('Enter the height of cylindrical part of 
acetylene bottle in meters '); 


Vol_top=(2*pi*r43)/3; % Calculating the 
volume of hemispherical top [m3] 
Vol_cyl=pi*r42*h; % Calculating the 


volume of cylindrical bottom [m3] 
Vol_total=Vol_top+Vol_cyl; % Calculating the 
total volume of acetylene bottle [m3] 

disp(' ') % Display blank line 
str = ['The volume of the acetylene bottle is ', 
num2str(Vol_total), ' cubic meters.']; 

disp(str); 

fo = fopen('output.txt', ‘w'); 

fprintf(fo,'The radius of acetylene bottle: %g 


meters \n', r); 

fprintf(fo, 'The height of cylindrical part of 
acetylene bottle: %g meters \n', h); 
fprintf(fo,'The volume of the acetylene bottle: %g 
cubic meters. \n', Vol_total); 

fclose(fo); 


Upon running the file, the output. txt file will display the following: 


The radius of acetylene bottle: 0.3 meters 

The height of cylindrical part of acetylene 
bottle: 1.5 meters 

The volume of the acetylene bottle: 0.480664 cubic 
meters. 


Loops 


In programming, a loop executes a set of code a specified number of times 
or until a condition is met. 


For Loop 


This loop iterates an index variable from an initial value using a specified 
increment to a final value and runs a set of code. The for loop syntax is the 
following: for loop_index=vector_statement code 

code end 


Example: 

Calculate y = cos(x) for —m < x < musing an increment of 7. for 
X=-p1:pi/4:pi y=cos(x); fprintf('%8.3f %8.2f 
Nit 5G, 2), 2enGd 


In the brief script above, x is the loop index that is initiated from —7 and 
incremented with 7 to a final value of 7. At the end of each increment, 

y = cos(z) is calculated and displayed with the fprintf command. This 
process continues until x = 7. 

From a previous exercise we know \n creates a new line when included in 
the fprintf command. Here, we also use %8.3f to specify eight spaces and 
three decimal places for the first variable x. Likewise %8.2f specifies the 
formatting for the second variable y but in this case, y is displayed with 
two decimal places. The result is the following: -3.142 -1.00 
-22550 =O, 71-1571 0,00) -O.705 05/1 0000 1,00 
Oevea Oba leo 0700 22350-0273 42 4 00) 

We can improve our code by adding formatting lines as follows: clear ; 
elle: fpranthi( > x cos(x)\n') % title row fprintit 
---------------- \n') % title row for x=-pi:pi/4:p1 
% 
loop_index=inital_value:increment_value: final_valu 
e y=cos(x); % code to calculate cos(x) 

fDFINET( "48.37 48.27 \n' ,x,y); % code tO print the 
output to screen end 

Screen output: X COS(X) ---------------- -3.142 -1.00 
e500 —O0(1 ol. 571 07007-00785 04/10, 000) 1.700 
O27S5. 0.71 12571 0,00 2. 3567-0, 71 3.142.-1- 00 


While Loop 


Like the for loop, a while loop executes blocks of code over and over again 
however it runs as long as the test condition remains true. The syntax of a 
while loop iswhile test_condition code ... code end 


Example: 
Using a while loop, calculate y = cos(z) for —7 < # < musing an 
increment of 7. 


This time we need to initialize the x value outside the loop and then state 
the test condition in the first line of the while loop. We also need to create 
an increment statement within the while loop: x=-p1; while x<=p1i 
V=COS(X) > TOriInth( 28-3 48.2 \n ,xX,V); x = x 4 
(pi/4); end 

The result is the same as that of the previous example: -3.142 -1.00 
220090 9077/1) S157 0,00 02 785 04/1) 0000 2.00 
0.785 0271) 1:574 0.002.356 —0.71-3,142 -1,00 

Now we can improve the code by adding extra formatting lines and 
comments: clear; clc; fprintf(' x cos(x)\n') % title 
row fprintf(' ---------------- \n') % title row x=- 
pl; % initiating the x value while x<=pi % stating 
the test condition y=cos(x); % calculating the 
Value of yy forantt(" 48.3f 48.2F Xn", x,y); % 
printing a and y xX = x + (pi/4); % iterating to 
the next step end 

The result should look the same as before. X COS(X) ------------ 
---- -3,142 -1.00 -2.356 -0.71 -1.571 0.00 -0.785 
O27. 0-000) 100707785. 0.740) 571 000 2. 356 7-0-7 1 
3.142 -1.00 


The 
diary 
Function 


Instead of writing a script from scratch, we sometimes solve problems in 
the Command Window as if we are using a scientific calculator. The steps 
we perform in this fashion can be used to create an m-file. For example, the 
diary function allows us to record a MATLAB session in a file and 
retrieve it for review. Reviewing the file and by copying relevant parts of it 
and pasting them in to an m-file, a script can be written easily. 


Typing diary at the MATLAB prompt toggles the diary mode on and off. 
As soon as the diary mode is turned on, a file called diary is created in the 
current directory. If you like to save that file with a specific name, say for 
example problem16, type 

>> diary problemié6.txt. 

A file named problem16.txt will be created. The following is the content of 
a diary file called problem16.txt. Notice that in that session, the user is 
executing the four files we created earlier. The user's keyboard input and the 
resulting display output is recorded in the file. The session is ended by 
typing diary which is printed in the last line. This might be useful to 
create a record of your work to hand in with a lab or to create the 
beginnings of an m-file. 


AcetyleneBottle 
Vol_total = 
0.4807 


AcetyleneBottleInteractive 

Enter the radius of acetylene bottle in meters .3 
Enter the height of cylinderical part of acetylene 
bottle in meters 1.5 


Vol_total = 
0.4807 


AcetyleneBottleInteractiveDisp 

This script computes the volume of an acetylene 
bottle 

Enter the radius of acetylene bottle in meters .5 
Enter the height of cylinderical part of acetylene 
bottle in meters 1.6 


The volume of the acetylene bottle is 
1.5184 


AcetyleneBottleInteractiveDisp1 
This script computes the volume of an acetylene 
bottle: 


Enter the radius of acetylene bottle in meters .9 
Enter the height of cylinderical part of acetylene 
bottle in meters 1.9 


The volume of the acetylene bottle is 6.3617 cubic 
meters. 
diary 


Style Guidelines 
Try to apply the following guidelines when writing your scripts: 


e Share your code or programs with others, consider adopting one of 
Creative Commons or GNU General Public License schemes 

e Include your name and contact info in the opening lines 

e Use comments liberally 

e Group your code and use proper indentation 

e Use white space liberally 

e Use descriptive names for your variables 

e Use descriptive names for your m-files 


Summary of Key Points 


1. A script is a file containing a sequence of MATLAB statements. Script 
files have a filename extension of .m. 

2. Functions such as input, disp and num2str can be used to make 
scripts interactive, 

3. Fopen, fprintf and fclose functions are used to create output 
files, 

4. A for loop is used to repeat a specific block of code a definite 
number of times. 


5. Awhile loop is used to repeat a specific block of code an indefinite 
number of times, until a condition is met. 

6. The diary function is useful to record a MATLAB command 
window session from which an m-file can be easily created, 

7. Various style guidelines covered here help improve our code. 


Problem Set 
Problem Set for Introductory Programming 
Exercise: 


Problem: 


Write a script that will ask for pressure value in psi and display the 
equivalent pressure in kPa with a statement, such as "The converted 
pressure is: ..." 


Solution: 


% This script converts pressures from psi to 
kPa % User is prompted to enter pressure in psi 
cle % Clear screen disp('This script converts 
pressures from psi to kPa:') disp(' ') % 
Display blank line psi=input('What is the 
pressure value in psi? '); kPa=psi*6.894757; % 
Calculating pressure in kPa disp(' ') % Display 
blank line str = ['The converted pressure is: 

+, MUMZSTY (KPa), ~ KPae |= Gisp( str); The script 
output is as follows: This script converts pressures 
from psi to kPa: What is the pressure value in 
psi? 150 The converted pressure is: 1034.2135 
kPa. 


Exercise: 


Problem: 


Write a script that generates a table of conversions from Fahrenheit to 
Celsius temperatures for a range and increment entered by the user, 
such as 


Enter the beginning temperature in F: 
Enter the ending temperature in F: 
Enter the increment value: 


Test your script with 20 the beginning Fahrenheit value, 200 the 
ending Fahrenheit value and 20 the increment. 


Solution: 


% This script generates a table of conversions 
% From Fahrenheit to Celsius temperatures clc % 
Clear screen disp('This script generates a 
table of conversions from Fahrenheit to 
Celsius’) disp(' ') % Display blank line 
lowerF=input('Enter the beginning temperature 
in F: '); upperF=input('Enter the ending 
temperature in F: '); increment=input('Enter 
the increment value: '); Fahrenheit= 

[ LowerF:increment:upperF]; % Creating a row 
vector with F values Celsius=5/9* (Fahrenheit - 
eZ) Yo CONVEFrEIng from F to 1G disp. ~ \os6 
Display blank line str = ['Fahrenheit Celsius 
'];% Displaying table header disp(str); % 
Tabulating results in two columns, ' is being 
used to transpose row to column 
disp([Fahrenheit' Celsius']) The script output is as 
follows: This script generates a table of 
conversions from Fahrenheit to Celsius Enter 
the beginning temperature in F: 20 Enter the 
ending temperature in F: 200 Enter the 


increment value: 20 Fahrenheit Celsius 20.0000 
-6.6667 40.0000 4.4444 60.0000 15.5556 80.0000 
26.6667 100.0000 37.7778 120.0000 48.8889 
140.0000 60.0000 160.0000 71.1111 180.0000 
82.2222 200.0000 93.3333 


Exercise: 
Problem: 
Pascal's Law states that pressure is transmitted undiminished in all 
directions throughout a fluid at rest. (See the illustration below). An 
initial force of 150 N is transmitted from a piston of 25 mm/2 to a 
piston of 100 mm/2. This force is progressively increased up to 200 N. 


Write a script that computes the corresponding load carried by the 
larger piston and tabulate your results. 


F4 my) 


P4 P95 


A simple hydraulic system. 


Solution: 


% This script computes the load carried by the 
larger piston in a hydraulic system clc % Clear 
screen disp('This script computes the load 
carried by the larger piston in a hydraulic 
system') disp(' ') % Display blank line 
initialF=150; finalF=200; increment=10; 
areail=25; area2=100; F1= 
[initialF:increment:finalF]; % Creating a row 
vector with Fi values F2=Fi*area2/areai; % 
Calculating F2 values disp(' ') % Display blank 
line str = | sFil F2 "| % Displayang, taple 
header disp(str); disp([F1' F2']) % Tabulating 
results in two columns, ' is being used to 
transpose row to column The script output is as follows: 
This script computes the load carried by the 
larger piston in a hydraulic system Fi F2 150 
600 160 640 170 680 180 720 190 760 200 800 


Exercise: 


Problem: 


Modify your script in previous problem so that the user provides the 
following input: 


Enter the initial force in N: 

Enter the final force in N: 

Enter the increment value: 

Enter the area of small piston in mm/2: 
Enter the area of big piston in mm/2: 


Test your script with 150, 200, 10, 25 and 100 with respect to each 
input variable. 


Solution: 


% This script computes the load carried by the 
larger piston in a hydraulic system clc % Clear 


screen disp('This script computes the load 
carried by the larger piston in a hydraulic 
system') disp(' ') % Display blank line 
initialF=input( ‘Enter the initial force in N: 
'); fainalF=input('Enter the final force in N: 
"); iIncrement=input('Enter the increment value: 
"); areai=input('Enter the area of small piston 
in mm42: '); area2=input('Enter the area of big 
piston in mm42: '); F1i= 

[initialF: increment: finalF]; % Creating a row 
vector with Fi values F2=Fi*area2/areali; % 
Calculating F2 values disp(' ') % Display blank 
line str = [' Fi F2 '|;% Displaying table 
header disp(str); disp([F1' F2']) % Tabulating 
results in two columns, ' is being used to 
transpose row to column The script output is as follows: 
This script computes the load carried by the 
larger piston in a hydraulic system Enter the 
initial force in N: 150 Enter the final force 
in N: 200 Enter the increment value: 10 Enter 
the area of small piston in mm42: 25 Enter the 
area of big piston in mm42: 100 F1 F2 150 600 
160 640 170 680 180 720 190 760 200 800 


Exercise: 


Problem: 
Write a script to solve the Stress-Strain problem in the Problem Set 
Solution: 


The m-file contains the following: % This script uses 
readings from a Tensile test and % Computes 
Strain and Stress values clc % Clear screen 
disp('This script uses readings from a Tensile 
test and') disp('Computes Strain and Stress 


values') disp(' ') % Display a blank line 
Specimen_dia=12.7; % Specimen diameter in mm % 
Load in kN Load_kN= 

O74 BOO. 779 A G7 ORO! 24045. 

27 662729 .397 32.68) 33,95; 34558 735.2252 % 
35./2;40..54;48..59°59.03°765.87'69..42>.... 
69.67;68.15;60.81]; % Gage length in mm 

length omm=(50 8; 50-8102, 50.8203 50.8305" 3. 
506406750. 8508; 50.8610); 50.8710 — 

5059016 75029270750. 9524-50. 97737 4.. 

510032; 51.8167 53.340" 55, 8380758. 420) 7 5. 
60..96761.468763.5;66.04]> % Calculate x- 
sectional area im m2 Cross_sectional_Area=pi/4* 
((Specimen_dia/1000)42); % Calculate change in 
length, initial lenght is 50.8 mm 
Delta_L=Length_mm-50.8; % Calculate Stress in 
MPa Sigma= 
(Load_kN./Cross_sectional_Area)*104(-3); % 
Calculate Strain in mm/mm 
Epsilon=Delta_L./50.8; str = ['Specimen 
diameter is ', num2str(Specimen_dia), ' mm.']; 
disp(str); Results=[Load_kN Length_mm Delta_L 
Sigma Epsilon]; % Tabulated results disp(' Load 
Length Delta L Stress Strain') disp(Results) 
After executed, the command window output is: This script 
uses readings from a Tensile test and Computes 
Strain and Stress values Specimen diameter is 
12.7 mm. Load Length Delta L Stress Strain 0 
50.8000 0 0 © 4.8900 50,8102 0.0102 38.6022 
0.0002 9.7790 50.8203 0.0203 77.1964 0.0004 
1456700 5058205 070305 115.8065 08-0006 19-5600 
50.8406 0.0406 154.4086 0.0008 24.4500 50.8508 
0:0508 19370108 0.0010 27.6200 50.8610 0.0610 
2157035190. 0012 29-3900 50-8710 0-07 11123270076 
80014 32.6800 5079016 OF 1016 257.9792 0.0020 
8329500) 5029270 021270 268-0047 ©0025 24-5800 


50.9524 071524 272.9780 0.0830 35.2200 50.9776 
Orly 732/78 .0302 070035 35.7200 Sil0032 0.2032 
281.9773 0.0040 40.5400 51.8160 1.0160 320.0269 
0.0200 48.3900 53.3400 2.5400 381.9955 0.0500 
29,0300 55.8800 5.0800 465.9388 071000 65.8700 
98.4200 7.6200 519.9844 0.1500 69.4200 60.9600 
10.1600 548.0085 0.2000 69.6700 61.4680 10.6680 
549, 962050. 2100687, 1500 63,5000 12; 7000 
537.9830 0.2500 60.8100 66.0400 15.2400 
480.0403 0.3000 


Exercise: 


Problem: 


Modify the script, you wrote above and plot an annotated Stress-Strain 
graph. 


Solution: 


Edited script contains the plot commands: % This script uses 
readings from a Tensile test and % Computes 
Strain and Stress values clc % Clear screen 
disp('This script uses readings from a Tensile 
test and') disp('Computes Strain and Stress 
values') disp(' ') % Display a blank line 
Specimen_dia=12.7; % Specimen diameter in mm % 
Load in kN Load_kN= 

KO 7 AO 97 (9 lA oy ORG? 244 oe oe. 

27.62; 29 39732-68733 ,95; 34568; 35.22) 2. 
35../2;40).54°48.359;°59,03°565..87'69..42>.... 
69.67;68.15;60.81]; % Gage length in mm 

Length omm=(50 28) 50281027 50.8703 50.8305. =... 
5078406; 50.8508; 50.8610; 50. 871i= =. 
5029016750. 927075029524 ° 50297 737.0. 

5100327 517816; 53.340" 55,880;58-420) 2. 

60.96: 61 .468763-5'66.04]; % Calculate x- 


sectional area im m2 Cross_sectional_Area=pi/4* 
((Specimen_dia/1000)42); % Calculate change in 
length, initial lenght is 50.8 mm 
Delta_L=Length_mm-50.8; % Calculate Stress in 
MPa Sigma= 
(Load_kN./Cross_sectional_Area)*104(-3); % 
Calculate Strain in mm/mm 
Epsilon=Delta_L./50.8; str = [ ‘Specimen 
diameter is ', num2str(Specimen_dia), ' mm.']; 
disp(str); Results=[Load_kN Length_mm Delta_L 
Sigma Epsilon]; % Tabulated results disp(' Load 
Length Delta L Stress Strain') disp(Results) % 
Plot Stress versus Strain plot(Epsilon, Sigma) 
title('Stress versus Strain Curve' ) 
Xlabel('Strain [mm/mm]') ylabel('Stress [mPa]') 
grid In addition to Command Window output, the following plot is 
generated: 


Stress versus Strain Curve 
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Exercise: 


Problem: 


Repeat Problem 2, this time using a combination of disp, fprintf 
commands and a for loop. 


Solution: 


The re-worked solution: % This script generates a table 
of conversions % From Fahrenheit to Celsius 
temperatures clear % removes all variables from 
the current workspace, clc % clears all input 
and output from the Command Window display, 
disp('This script generates a table of 
conversions from Fahrenheit to Celsius') disp(' 
") % Display blank line lowerF=input('Enter the 
initial temperature in F: '); 
upperF=input('Enter the final temperature in F: 
"); LIncrement=input('Enter the increment value: 
‘yo disp(* *) % Display blank line 

fprintf( ‘Fahrenheit Celsius\n') % title row 
fprintf('------------------ \n') % title row for 
Fahrenheit=[lowerF:increment:upperF]; % 
Creating a row vector with F values 
Celsius=5/9* (Fahrenheit-32); % Converting from 
FLO, C fprinkn( Yesah %e8.o0 

\n', Fahrenheit,Celsius); % Tabulating results 
in two columns end 

After executed, the command window output is: This script 
generates a table of conversions from 
Fahrenheit to Celsius Enter the initial 
temperature in F: 20 Enter the final 
temperature in F: 200 Enter the increment 
value: 20 Fahrenheit Celsius ------------------ 
20.000 -6.667 40.000 4.444 60.000 15.556 80.000 


26.667 100.000 37.778 120.000 48.889 140.000 
60,000, 160 000 7 tii 180) 000 622222 200,000 
SiS) rere! 


Exercise: 


Problem: Repeat Problem 7, this time using a while loop. 
Solution: 


The re-worked solution: % This script generates a table 
of conversions % From Fahrenheit to Celsius 
temperatures clear % removes all variables from 
the current workspace, clc % clears all input 
and output from the Command Window display, 
disp('This script generates a table of 
conversions from Fahrenheit to Celsius') disp(' 
") % Display blank line lowerF=input('Enter the 
initial temperature in F: '); 
upperF=input('Enter the final temperature in F: 
"); increment=input('Enter the increment value: 
"); disp(' ') % Display blank line 

fprintf( ‘Fahrenheit Celsius\n') % title row 
fprintf('------------------ \n') % title row 
Fahrenheit=lowerF; while Fahrenheit<=upperF 
Celsius=5/9* (Fahrenheit-32); % Converting from 
E COs hOeiinen( woot, oo 231i 

\n', Fahrenheit,Celsius); % Tabulating results 
in two columns Fahrenheit=Fahrenheit+increment; 
end 

After executed, the command window output is: This script 
generates a table of conversions from 
Fahrenheit to Celsius Enter the initial 
temperature in F: 20 Enter the final 
temperature in F: 200 Enter the increment 
value: 20 Fahrenheit Celsius ------------------ 


20.000 -6.667 40.000 4.444 60.000 15.556 80.000 
26,6607 100,000 37.778 120,000 48.839 140,000 
60,000) 160,000 71,111 1807 000182.222 200,000 
93.333 


Interpolation 
Interpolation with MATLAB. 


linear T mterpélaction 


www, hetemeal.com 


Linear interpolation is one of the most common techniques for estimating 
values between two given data points. For example, when using steam 
tables, we often have to carry out interpolations. With this technique, we 
assume that the function between the two points is linear. MATLAB has a 
built-in interpolation function. 


The 
interp1 


Function 


Give an x-y table, y_new = interp1(x, y, X_new) interpolates to 
find y_new. Consider the following examples: 


Example: 

To demonstrate how the interp1 function works, let us create an x-y 
table with the following commands; xX = 0:5; y = 

[0,20,60,68, 77,110]; To tabulate the output, we can use 
[x',¥  ] Wie@resniiisans = 0 0 1 20 2 60 3 68 4 77 5 
110 Suppose we want to find the corresponding value for 1.5 or 
interpolate for 1.5. Using y_new = interp1(x,Y,X_new) syntax, 
let us type in: y_new=interp1i(x,y,1.5) y_new = 40 


Example: 

The table we created above has only 6 elements in it and suppose we need 
a more detailed table. In order to do that, instead of a single new x value, 
we can define an array of new x values, the interp1 function returns an 
array of new y values: new_X = 0:0.2:5; new_y = 

interp1(x, y,new_x); Let's display this table [New_x',new_y' | 
The result is ans = 0 0 0.2000 4.0000 0.4000 8.0000 
-6000 12.0000 0:S000 16.0000 1.0000 20.0000 

.2000 28.0000 1.4000 36.0000 1.6000 44.0000 
70000) 52,0000 2.0080 60,0000) 2.2000 61,6000 
.4000 63.2000 2.6000 64.8000 2.8000 66.4000 
.0000 68.0000 3.2000 69.8000 3.4000 71.6000 
.6000 73.4000 3.8000 75.2000 4.0000 77.0000 
.2000 83.6000 4.4000 90.2000 4.6000 96.8000 
-8000 103.4000 5.0000 110-0000 


BRBRWWNHERE O 
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Example: 
Using the table below, find the internal energy of steam at 215 °C and the 
temperature if the internal energy is 2600 kJ/kg (use linear interpolation). 


Temperature [C] Internal Energy [kJ/kg] 


100 2506.7 
150 2582.8 
200 2658.1 
250 2733.7 
300 2810.4 
400 2967.9 
500 3131.6 


An extract from Steam Tables 


First let us enter the temperature and energy values temperature = 
[100, 150, 200, 250, 300, 400, 500]; energy = 
[250637,. 256200, 2oocu 2iocur-. 2000.4 2967.9) 
3131.6]; [temperature',energy'] returmnsans = 
1e0e700S * 0.1000) 2.5067 0.1500) 2753822 0.2000 
Z.0501 0, 2500 2773427 ©,3008 2.8104 0.4000) 2.9679 
0.5000 3.1316 issue the following for the first question: 
new_energy = interpi(temperature, energy, 215) returns 
new_energy 2.6808e+003 Now, type in the following for the 
second question: new_temperature = 

interpi(energy, temperature, 2600) returns 
new_temperature = 161.4210 


Summary of Key Points 


1. Linear interpolation is a technique for estimating values between two 
given data points, 

2. Problems involving steam tables often require interpolated data, 

3. MATLAB has a built-in interpolation function. 


Problem Set 
Problem Set for Interpolation with MATLAB 
Exercise: 


Problem: 
Determine the saturation temperature, specific liquid enthalpy, specific 


enthalpy of evaporation and specific enthalpy of dry steam at a 
pressure of 2.04 MPa. 


Saturation 
Pressure ‘Temperature hy heg hg 
[MN/m7] [Cc] [kJ/kg] [kJ/kg] [kJ/kg] 
251 214.9 920.0 1878.2 2798.2 
2.0 212.4 908.6 1888.6 2797.2 


An extract from steam tables 


Solution: 


MATLAB solution is as follows; >> pressure=[2.1 2.0]; >> 
sat_temp=[214.9 212.4]; >> h_f=[920 908.6]; >> 
hofo=(1878.2 1863.6], == Neg=(279e8.2 2797.2]; 
>> sat_temp_new=interp1(pressure, sat_temp, 2.04) 
sat_temp_new = 213.4000 >> 
h_f_new=interpi(pressure,h_f,2.04) h_f_new = 
913.1600 >> 
h_fg_new=interpi(pressure,h_fg,2.04) h_fg_new = 
1.8844e+003 >> 


h_g_new=interp1(pressure,h_g,2.04) h_g_new = 
227976270038 


Exercise: 


Problem: 


The following table gives data for the specific heat as it changes with 
temperature for a perfect gas. (Data available for download). [footnote] 


Temperature [F] Specific Heat [BTU/AbmF] 
25 0.118 
50 0.120 
75 0.123 
100 0.125 
125 0.128 
150 0.131 


Change of specific heat with temperature 


Using interp1 function calculate the specific heat for 30 F, 70 F and 
145 F. 

Thermodynamics and Heat Power by Kurt C. Rolle, Pearson Prentice 
Hall. © 2005, (p.19) 


Solution: 


MATLAB solution is as follows: >> temperature= 
(255503757 10071257150] temperature = 25°50 75 
100 125 150 >> specific_heat= 

jedte sl 20) 123) .025° 128° 131 specie heat = 
Orts0) O,i700 07 1220)0. 1250) Oni 80 O.a3s10) Se 
specific_heatAt30=interp1(temperature, specific_ 
heat,30) specific_heatAt30 = 0.1184 >> 
specific_heatAt70=interp1( temperature, specific_ 
heat, 70) specific_heatAt70 = 0.1224 >> 
specific_heatAt145=interp1i(temperature, specific 
_heat,145) specific_heatAt145 = 0.1304 


Exercise: 


Problem: 


For the problem above, create a more detailed table in which 
temperature varies between 25 and 150 with 5 F increments and 
corresponding specific heat values. 


Solution: 


MATLAB solution is as follows: >> 
neworemperatunpe=25:5:1507 >> 
new_specific_heat=interp1i(temperature, specific_ 
heat,new_temperature); >> 
[new_temperature’',new_specific_heat'] a 
25.0000 6.1180 3070000 0.1184 35.0000 6 
40.0000 0.1192 45.0000 0.1196 50.0000 0 
55,0000) 0 -1206 60,0000 0.1212 65,0000 6 
70.0000 ©.1224 75-0000 ©.-1230 80.0000 0 
85,0080 G.1238 90,0000) 0, 1242 95.0000) 6.1246 
HOO 0G00 01250 2105-00000 1256. 110.0000 0: 1262 
d15.00800) 071268 120.0000 0.11274) 125.0000 0.1280 
130.0000 0.1286 135.0000 021292 14070000 0.1298 
145.0000 0) 1364 150.0000 0.1316 


Za! 


Exercise: 


Problem: 


During a 12-hour shift a fuel tank has varying levels due to 
consumption and transfer pump automatically cutting in and out to 
maintain a safe fuel level. The following table of fuel tank level versus 
time (Data available for download) is missing readings for 5 and 9 
AM. Using linear interpolation, estimate the fuel level at those times. 


Time [hours, AM] Tank level [m] 
1:00 1.5 
2:00 Py, 
3:00 22 
4:00 29 
5:00 ? 
6:00 2.0 
7:00 20 
8:00 2:3 
9:00 ? 
10:00 2.0 


11:00 1.8 


Time [hours, AM] Tank level [m] 
12:00 1.3 


Fuel tank level versus time 


Solution: 


SS eiie= (deen ou4 aon oO Lonel 

tank 1evel=[io5 2.7 273° 2.9 2,6 2,5, 253 2.0 1.8 
dol 

tank_level_at_5=interpi(time, tank_level,5) 
tank_level_at_5 = 2.7500 >> 
tank_level_at_9=interpi(time, tank_level,9) 
tank_level_at_9 = 2.1500 


Computing the Area Under a Curve 
Basic Numerical Integration with MATLAB 


Area under a Curve 
T rapezoidal Rule 
Arorymoas Fone See. 
q Integra) Func +ian 
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This chapter essentially deals with the problem of computing the area under 
a curve. First, we will employ a basic approach and form trapezoids under a 
curve. From these trapezoids, we can calculate the total area under a given 
curve. This method can be tedious and is prone to errors, so in the second 
half of the chapter, we will utilize a built-in MATLAB function to carry out 
numerical integration. 


A Basic Approach 


There are various methods to calculating the area under a curve, for 
example, Rectangle Method, Trapezoidal Rule and Simpson's Rule. The 
following procedure is a simplified method. 


Consider the curve below: 


At Av Ar At At 


Numerical 
integration 


Each segment under the curve can be calculated as follows: 
Equation: 


1 1 1 
2 (yo + y1)Ax + m1 (yr + y2)Ax + o (yo + y3)Ax 


Therefore, if we take the sum of the area of each trapezoid, given the limits, 
we calculate the total area under a curve. Consider the following example. 


Example: 


Given the following data, plot an x-y graph and determine the area under a 
curve between x=3 and x=30 


Index x [m] y [IN] 


Index x [m] y [IN] 


iE 3 27.00 

Z 10 14.50 

3 15 9.40 

4 20 6.70 

5 25 5.30 

6 30 4.50 
Data Set 


First, let us enter the data set. For x, issue the following command x= 
[3,10,15, 20, 25, 30] ; Raationay= 
[27,14.5,9.4,6.7,5.3,4.5] ; Mime [ xX", y' | Byitesnlll 
see the following tabulated result. Here we transpose row vectors with ' 
and displaying them as columns: 


ans = 
3.0000 27.0000 
10.0000 14.5000 
15.0000 9.4000 
20.0000 6. 7000 
25.0000 5.3000 
30.0000 4.5000 


Compare the data set above with the given information in the question. 
To plot the data type the following: 


plot(x,y),title( 'Distance-Force 
Graph'),xlabel('Distance[m]'),ylabel('Force[N]'),g 
rid 


The following figure is generated: 
Distance-Force Graph 


15 


Forea(N] 


Distance-Force Graph 


To compute dx for consecutive x values, we will use the index for each x 
value, see the given data in the question.: 


dx=[x(2)-x(1),x(3)-x(2),x(4)-x(3),x(5)-x(4),x(6) - 
x(5)]; 


dy is computed by the following command: 


dy=[0.5*(y(2)+y(1)),0.5*(y(3)ty(2)),0.5* 
(y(4)+y(3)),0.5*(y(5)t+y(4)),0.5*(y(6)+y(5) ) J; 


dx and dy can be displayed with the following command: [dx',dy' ]. 
The result will look like this: 


[dx',dy' ] 
ans = 


. 0000 20.7500 
. 0000 11.9500 
. 0000 8.0500 
. 0000 6.0000 
. 0000 4.9000 


o1 O1 O1 O1 N 


Our results so far are shown below 


x [m] y [IN] dx [m] dy [N] 
3 27.00 

10 14.50 7.00 20.75 
15 9.40 5.00 11.95 
20 6.70 5.00 8.05 
25 5.30 5.00 6.00 
30 4.50 5.00 4.90 


x, y and corresponding differential elements 


If we multiply dx by dy, we find da for each element under the curve. The 
differential area da=dx*dy, can be computed using the 'term by term 
multiplication’ technique in MATLAB as follows: 


da=dx. *dy 
da = 

145.2500 59.7500 40.2500 30.0000 24.5000 
Each value above represents an element under the curve or the area of 
trapezoid. By taking the sum of array elements, we find the total area under 
the curve. 
sum(da) 
ans = 


299.7500 


The following illustrates all the steps and results of our MATLAB 
computation. 


x [m] y [N] dx [m] dy [N] dA [Nm] 
3 27.00 

10 14.50 7.00 20.75 145.25 
15 9.40 5.00 11.95 59.75 

20 6.70 5.00 8.05 40.25 

25 5.30 5.00 6.00 30.00 


30 4.50 5.00 4.90 24.50 


x [m] y [IN] dx [m] dy [N] dA [Nm] 
2) 3 


Computation of the approximate area under a curve 


The Trapezoidal Rule 


Sometimes it is rather convenient to use a numerical approach to solve a 
definite integral. The trapezoid rule allows us to approximate a definite 
integral using trapezoids. 


The 
trapz 
Command 


Z = trapz(Y) computes an approximation of the integral of Y using the 
trapezoidal method. 


Now, let us see a typical problem. 


Example: 


Given Area = HE : az” d x, an analytical solution would produce 39. Use 
trapz command and solve it 


1. Initialize variable x as a row vector, from 2 with increments of 0.1 to 
max=2:.1:5; 

2. Declare variable y as y=x2;. Note the following error prompt: ??? 
Error using ==> mpower Inputs must be a scalar 


and a square matrix. This is because x is a vector quantity 
and MATLAB is expecting a scalar input for y. Because of that, we 
need to compute y as a vector and to do that we will use the dot 
operator as follows: y=x .42;. This tells MATLAB to create vector y 
by taking each x value and raising its power to 2. 

3. Now we can issue the following command to calculate the first area, 
the output will be as follows: 


areail=trapz(x,y) 
areal = 
39.0050 


Notice that this numerical value is slightly off. So let us increase the 
number of increments and calculate the area again: 


K=24 0075) 
y=x.A2; 
area2=trapz(x,y) 
area2 = 
39.0001 
Yet another increase in the number of increments: 
X=2:.001:57 
y=x.A2; 
area3=trapz(x,y) 


area3 = 


39.0000 


Example: 
Determine the value of the following integral: 
fo sin(z) da 


1. Initialize variable x as a row vector, from 0 with increments of pi/100 
to pi: Xx=0:pi1/100:pi; 

2. Declare variable y as y=Sin(x); 

3. Issue the following command to calculate the first area, the output 
will be as follows: 


areail=trapz(x,y) 
areal = 
1.9998 


let us increase the increments as above: 


X=0:p1i/1000: pi; 
y=sin(x); 
area2=trapz(x,y) 


area2 = 


2.0000 


Example: 

A gas expands according to the law, EN =o Initially, the pressure is 100 
kPa when the volume is 1 m®. Write a script to compute the work done by 
the gas in expanding to three times its original volume[footnote]. 

O. N. Mathematics: 2 by J. Dobinson, Penguin Library of Technology. © 
1969, (p. 184) 

Recall that PV diagrams can be used to estimate the net work performed by 
a thermodynamic cycle, see Wikipedia or we can use definite integral to 


compute the work done (WD) as follows: 
Equation: 


wD = [pds 


If we rearrange the expression pressure as a function of volume, we get: 
Equation: 


By considering the initial state, we can determine the value of c: 
Equation: 


ec = 100~x1!4 


100 
From the equation and the equation above, we can write: 
Equation: 
100 
P — 
yl4 


By inserting P in WD, we get: 
Equation: 


For MATLAB solution, we will consider P_as a function of V and WD. 
Now, let us apply the three-step approach we have used earlier: 


1. Initialize variable volume as a row vector, from 1 with increments of 
OU0MeB32V=1:0.001:3; 

2. Declare variable pressure as p=100./v.41.4; 

3. Use the trapz function to calculate the work done, the output will be 
as follows: 


WorkDone=trapz(v,p) 
WorkDone = 
88.9015 


These steps can be combined in an m-file as follows: 


cle 
disp('A gas expands according to the law, 
pvA1.4=C') 


disp('Initial pressure is 100 kPa when the volume 
is 1 m3') 

disp('Compute the work done by the gas in 
expanding’ ) 

disp('To three times its original volume' ) 


disp(' ') % Display blank line 
v=1: .001:3; % Creating a row vector for 
volume, v 


p=100./(v.41.4); % Computing pressure for 
volume 

WorkDone=trapz(v,p) % Integrating p*dv over 1 to 3 
Cubic meters 


Example: 
A body moves from rest under the action of a direct force given by 


ie — where x is the distance in meters from the starting point. Write a 


script to compute the total work done in moving a distance 10 m.[footnote | 
O. N. Mathematics: 2 by J. Dobinson, Penguin Library of Technology. © 
1969, (p. 183) 

Recall that the general definition of mechanical work is given by the 
following integral, see Wikipedia: 

Equation: 


wo = { Fdz 


Therefore we can write: 
Equation: 


10 
15 
wp = | dz 
0 zr+3 


Applying the steps we followed in the previous examples, we write: 


cle 

disp('A body moves from rest under the action of a 
direct force given' ) 

disp('by F=15/(x+3) where x is the distance in 
meters' ) 

disp('From the starting point. ') 

disp('Compute the total work done in moving a 
distance 10 m.') 

disp(' ') % Display blank line 
xX=0: .001:10; % Creating a row vector for 
distance, xX 

F=15./(x+3); % Computing Force for x 
WorkDone=trapz(x,F) % Integrating F*dx over 0 to 
10 meters. 


The output of the above code is: 


A body moves from rest under the action of a 
direct force given 

by F=15/(x+3) where x is the distance in meters 
From the starting point. 

Compute the total work done in moving a distance 
10 Mm. 


WorkDone = 


Zi oo ou 


The 
integral 


Function 


As we have seen earlier, tf apZ gives a good approximation for definite 
integrals. The integral function streamlines numerical integration even 
further. Before we learn about integral function, first we will look at 
anonymous functions. 


Anonymous Functions 


An anonymous function is a function that can be defined in the command 
window (i.e. it does not need to be stored in a program file). Anonymous 
functions can accept inputs and return outputs, just as standard functions do 
such as sqrt(X) or log(X). 


To define an anonymous function, first we create a handle with @( x) and 
type in the function: myFfunction=@(x) x42+1. 


If you want to evaluate myfunction at 1, just type in 
a=myfunction(1) at the command window and you get the result of 2. 


Syntax for integral 

To evaluate an integral from a minimum to a maximum value, we specify a 
function and its minimum and maximum Z = 
integral( fun, xmin, xmax). 


Example: 
Given y=xX‘2, evaluate the integral from x=2 to x=5 as we have done it 
with trapz command. 


1. Define function myfunction=@(x) x.42; 

2. Apply the syntax to myfunction as follows Z = 
integral(myfunction, 2,5) 

3. You should get a result of Z = 39. 


Note: Notice that, unlike in 
trapz 


example, we did not need to define a vector and change the increments to 
get an accurate result. 


Summary of Key Points 


1. In its simplest form, numerical integration involves calculating the 
areas of segments that make up the area under a curve, 

2. MATLAB has built-in functions to perform numerical integration, 

3.Z = trapz(Y) computes an approximation of the integral of Y 
using the trapezoidal method. 

4. Anonymous functions are inline statements that we can define with 
@(x), 

5.Z = integral( fun, xmin, xmax) numerically integrates 
function fun from Xmin to xmax. 


Problem Set 
Problem Set for Numerical Integration with MATLAB 
Exercise: 


Problem: 


Let the function y defined by y = cos(z). Plot this function over the 
interval [-pi,pi]. Use numerical integration techniques to estimate the 
integral of y over [0, pi/2] using increments of pi/10 and pi/1000. 


Solution: 


1. Plotting: x=-pi:pi/100:p1i; y=cos(x); 
plot Cc y), title Graph or 
y=cos(xX)'),xlabel('x'),ylabel('y'), grid 

2. Area calculation 1: >> x=0:pi/10:pi/2; >> y=cos(x); 
>> areai=trapz(x,y) areal = 0.9918 

3. Area calculation 2: >> x=0:p1i/1000:p1i/2; >> 
y=cos(X); >> area2=trapz(x,y) area2 = 1.0000 


Exercise: 
Problem: 
Let the function y defined by y = 0.04x? — 2.132 + 32.58. Plot this 
function over the interval [3,30]. Use numerical integration techniques 
to estimate the integral of y over [3,30]. 


Solution: 


1. Plotting:>> x=3:.1:30; >> y=0.04* 
(X28 2 = 2513." xX+32,58) >> ploux, yy), 


title( ‘Graph of ... y=.04* 
(X42)-2.13*x+32-58" ), xllabel('’x"), ylabel(C’y” ) 
,grid 


2. Area calculation: >> area=trapz(x,y) area = 
290.3868 


Exercise: 


Problem: 


A 2000-liter tank is full of lube oil. It is known that if lube oil is 
drained from the tank, the mass flow rate will decrease from the 
maximum when the tank level is at the highest. The following data 
were collected when the tank was drained. 


Time [min] Mass Flow [kg/min] 
0 50.00 

5 48.25 

10 46.00 

15 42.50 

20 37.50 

25 30.50 

30 19.00 

35 9.00 

Data 


Write a script to estimate the amount of oil drained in 35 minutes. 


Solution: 


clc t=linspace(0,35,8) % Data entry for time 
[min] m=[50 48.25 46 42.5 37.5 30.5 19 9] % 
Data entry for mass flow [kg/min] % Calculate 
time intervals dt=[t(2)-t(1),t(3)-t(2),t(4)- 
CE) pee SSI Sy RS CTO Ie ye el te = 
t(7)] % Calculate mass out dm=[0.5* 
(m(2)+m(1)),0.5*(m(3)+m(2)),0.5* 
CNM) Os (NCS ote cI 4) nO or 
(m(6)+m(5)),0.5*(m(7)+m(6)),0.5*(m(8)+m(7))] % 
Calculate differential areas da=dt.*dm; % 
Tabulate time and mass flow [t',m'] % Tabulate 
time intervals, mass out and differential areas 
[dt',dm',da'] % Calculate the amount of oil 
drained [kg] in 35 minutes 01il_Drained=sum(da) 
The output is:ans = 0 50.0000 5.0000 48.2500 

16 @€080 46.0000) 15-0000) 42-5000 20,0000 37-5000 
25.0000 30.5000 30.0000 19.0000 35.0000 9.0000 
ans = 5.0000 49.1250 245.6250 5.0000 47.1250 
235.6250 5.0000) 44° 2500) 22172508 5.0000 40.0000 
200.0000 520000 24.0000 176.0008 5:-0000 24.7500 
123.7500 5.0000 14.0000 70.0000 011 Drained = 
i. 266367003 


Exercise: 


Problem: 


A gas is expanded in an engine cylinder, following the law PV!=c. 
The initial pressure is 2550 kPa and the final pressure is 210 kPa. If the 
volume at the end of expansion is 0.75 m?, compute the work done by 
the gas. [footnote] 

Applied Heat for Engineers by W. Embleton and L Jackson, Thomas 
Reed Publications. © 1999, (p. 80) 


Solution: 


cle disp('A gas is expanded in an engine 
cylinder, following the law PV‘1.3=c' ) 
disp('The initial pressure is 2550 kPa and the 
final pressure is 210 kPa.') disp('If the 
volume at the end of expansion is 0.75 m3,') 
disp('Compute the work done by the gas.') 
disp(' ') % Display blank line n=1.3; P_1=2550; 
% Initial pressure P_f=210; % Final pressure 
V_f=.75; % Final volume V_1i=(P_T* 
(V_f4n)/P_1)4(1/n); % Initial volume 
c=P_f*V_fA4n; v=V_1:.001:V_f; % Creating a row 
vector for volume, v p=c./(v.4n); % Computing 
pressure for volume WorkDone=trapz(v,p) % 
Integrating p*dv The outputis: A gas 1S expanded in 
an engine cylinder, following the law PV41.3=c 
The initial pressure is 2550 kPa and the final 
pressure is 210 kPa. If the volume at the end 
of expansion is 9.75 m3, Compute the work done 
by the gas. WorkDone = 409.0666 


Exercise: 


Problem: 


A force F acting on a body at a distance s from a fixed point is given 
by # = 3s + _ Write a script to compute the work done when the 
body moves from the position where s=1 to that where s=10. 
[footnote | 

O. N. Mathematics: 2 by J. Dobinson, Penguin Library of Technology. 
© 1969, (p. 213) 


Solution: 


cle disp('A force F acting on a body at a 
distance s from a fixed point is given by') 
disp('F=3*s+(1/(s42)) where s is the distance 
in meters') disp('Compute the total work done 


in moving') disp('From the position where s=1 
to that where s=10.') disp(' ') % Display blank 
line s=1:.001:10; % Creating a row vector for 
Gistance, Ss F=3.°sr(i2/(s.%2)); % Computing 
Force for s WorkDone=trapz(s,F) % Integrating 
F*ds over 1 to 10 meters. The output is: A force F 
acting on a body at a distance s from a fixed 
point is given by F=3*s+(1/(s42)) where s is 
the distance in meters Compute the total work 
done in moving From the position where s=1 to 
that where s=10. WorkDone = 149.4000 


Exercise: 


Problem: 


The pressure p and volume v of a given mass of gas are connected by 
the relation (p + =) (v — b) = k where a, b and k are constants. 
Express p in terms of v, and write a script to compute the work done 
by the gas in expanding from an initial volume to a final volume. 
[footnote | 


Test your solution with the following input: 

a: 0.01 

b: 0.001 

The initial pressure [kPa]: 100 

The initial volume [m3]: 1 

The final volume [m3]: 2 

O. N. Mathematics: 2 by J. Dobinson, Penguin Library of Technology. 
© 1969, (p. 212) 


Solution: 


clc % Clear screen disp('This script computes 
the work done by') disp('The gas in expanding 
from volume vi to v2') disp(' ') % Display 
blank line a=input( ‘Enter the constant a: '); 
b=input('Enter the constant b: '); 


p_i=input('Enter the initial pressure [kPa]: 

'); v_i=input('Enter the initial volume [m3]: 
'); v_f=input('Enter the final volume [m3]: '); 
k=(plit(a/ (v_1%2))*(v-1-b) ); % Calculating 
constant k v=v_1:.001:v_f; % Creating a row 
vector for volume p=(k./(v-b))-(a./(v.42)); % 
Computing pressure for volume 
WorkDone=trapz(v,p); % Integrating p*dv disp(' 
') % Display blank line str = ['The work done 
by the gas in expanding from‘, 

MUMNZSEP (VoL) ae. MNS EO NUNZStr (ven), 7 me 
1s ', num2str(WorkDone), * kW. "|; dasp(str); 
The output is: This script computes the work done by 
The gas in expanding from volume vi to v2 Enter 
the constant a: .01 Enter the constant b: .001 
Enter the initial pressure [kPa]: 100 Enter the 
initial volume [m3]: 1 Enter the final volume 
[m3]: 2 The work done by the gas in expanding 
from 1 m3 to 2 m2 1s 69).38667 kW. 


Regression Analysis 
Basic curve fitting and linear regression. 


line ar Regression 
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What is Regression Analysis? 


Suppose we calculate some variable of interest, y, as a function of some 
other variable x. We call y the dependent variable and x the independent 
variable. For example, consider the data set below, taken from a simple 
experiment involving a vehicle, its velocity versus time is tabulated. In this 
case, velocity is a function of time, thus velocity is the dependent variable 
and the time is the independent variable. 


Time [s] Velocity [m/s] 


Time [s] Velocity [m/s] 


0 20 
10 39 
20 67 
30 89 
40 111 
50 134 
60 164 
70 180 
80 200 


Vehicle velocity versus time. 


In its simplest form regression analysis involves fitting the best straight line 
relationship to explain how the variation in a dependent variable, y, depends 
on the variation in an independent variable, x. In our example above, once 
the relationship (in this case a linear relationship) has been estimated we 
can produce a linear equation in the following form: 

Equation: 


y=mrzet+n 


And once an analytic equation such as the one above has been determined, 
dependent variables at intermediate independent values can be computed. 


Performing Linear Regression 


Regression analysis with MATLAB is easy. The MATLAB Basic Fitting 
GUI allows us to interactively to do "curve fitting" which is a method to 
arrive at the best "straight line" fit for linear equations or the best curve fit 
for a polynomial up to the tenth degree. The procedure to perform a curve 


fitting with MATLAB is as follows: 


1. Input the variables, 


2. Plot the data, 


3. Initialize the Basic Fitting GUI, 


4. Select the desired regression analysis parameters. 


Example: 


Using the data set above, determine the relationship between velocity and 


time. 


First, let us input the variables (Workspace > New variable) as shown in 
the following figures. 


MATLAB 7.11.0 (R2010b) 
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A new variable is created in the Workspace. 
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Select a file to view details 


New variables are entered in the Variable Editor. 


Second, we will plot the data by typing in plot (time, velocity) at 
the MATLAB prompt. The following plot is generated, select Tools > 
Basic Fitting: 


I lolx] 


File Edit edie Insert | Tools Desktop Window Help 
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A plot is generated in Figure 1. The Basic Fitting tool can be 
initialized from Tools > Basic Fitting. 


In the "Basic Fitting" window, select "linear" and "Show equations". The 
best fitting linear line along with the corresponding equation are displayed 
on the plot: 
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Basic Fitting window is used to select the desired regression analysis 
parameters. 


Now let us do another curve fitting and obtain an equation for the function. 
Using that equation, we can evaluate the function at a desired value with 
polyval. 


Example: 

The following is a collection of data for an iron-constantan thermocouple 
(data available for download). [footnote] 

Engineering Fundamentals and Problem Solving by Arvid R. Eide, Roland 
Jenison, Larry L. Northup, Steven K. Mikelson , McGraw-Hill Higher 
Education. © 2007 p.114 


Temperature [C] 
50 
100 
150 
200 
300 
400 
500 
600 
700 
800 
900 
1000 


Temperature [C] vs Voltage [mV] 


Voltage [mV] 
276 
6.7 
8.8 
IGE 
17.0 
ae: 
26 
32.5 
a7 
41 
48 


Daud 


a. Plot a graph with Temperature as the independent variable. 
b. Determine the equation of the relationship using the Basic Fitting 


tools. 


c. Estimate the Voltage that corresponds to a Temperature of 650 C and 


MENG: 


We will input the variables first: 


Temp= 


[50; 100; 150; 200; 300; 400; 500; 600; 700; 800; 900; 1000 | 


Voltage= 


([2e67Ga7- 8.0 2 ai 22557 26752.57 377 | Ad AG 55.2] 


To plot the graph, type in: 


plot(Temp, Voltage) 


We can now use the Plot Tools and Basic Fitting settings and determine the 


equation: 


A) Figure 1 TE ee ites we cite we TEES basic Fitting. 
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Basic Fitting window is used to select the desired regression analysis 


parameters. 


By clicking the right arrow twice at the bottom right corner on the Basic 
Fitting window, we can evaluate the function at a desired value. See the 
figure below which illustrates this process for the temperature value 1150 
C. 
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Estimating the Voltage that corresponds to a Temperature of 1150 C. 


Now let us check our answer with a technique we learned earlier. As 
displayed on the plot, we have obtained the following equation: 

y = 0.052831z + 0.67202 This equation can be entered as polynomial 
and evaluated at 650 and 1150 as follows: 


>> p=[0.052831,0.67202] 
p= 


0.0528 0.6720 


>> polyval(p, 1150) 
ans = 


61.4277 


Summary of Key Points 


1. Linear regression involves fitting the best straight line relationship to 
explain how the variation in a dependent variable, y, depends on the 
variation in an independent variable, x, 

2. Basic Fitting GUI allows us to interactively perform curve fitting, 

3. Some of the plot fits available are linear, quadratic and cubic functions, 

4. Basic Fitting GUI can evaluate functions at given points. 


Problem Set 
Problem Set for Regression Analysis with MATLAB 
Exercise: 


Problem: 
Using the following experimental values [footnote], plot a distance- 


time graph and determine the equation, relating the distance and time 
for a moving object. 


Distance [m] Time [s] 
0 0 

24 5 

48 10 

72 15 

96 20 


Experimental data. 


Engineering Science by E. Hughes and C. Hughes, Longman © 1994, 
(p. 375) 


Solution: 


Data can be entered as follows: distance=[0 24 48 72 96]; 
time=[0 5 10 15 20]; wecannow plot the data by typing in 
plot(time, distance);title('Distance-Time 

Graph' );xlabel( 'time');ylabel('distance' ); atthe 


MATLAB prompt. The following plot is generated, select Tools > 
Basic Fitting: 
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As shown above, the relationship between distance and time is: 
y =4.82 —1.7x 10-4 

or 

Distance = 4.8 Time — 1.7 x 10-4 


Exercise: 


Problem: 


Using the data set below, determine the relationship between 
temperature and internal energy. 


Temperature [C] Internal Energy [kJ/kg] 


100 2506.7 


Temperature [C] Internal Energy [kJ/kg] 


150 2582.8 
200 2658.1 
250 273347 
300 2810.4 
400 2967.9 
500 3131.6 


An extract from Steam Tables 


Solution: 


Data can be entered as follows: temperature = [100, 150, 
200, 250, 300, 400, 500]; energy = [2506.7, 
2502.0, 2058.1, 2738.7, 200.4, 2967.9, 
3131.6]; we can now plot the data by typing in 
plot(temperature, energy);title('temperature vs. 
energy');xlabel('temperature' );ylabel('energy' ) 
; at the MATLAB prompt. The following plot is generated, select 
Tools > Basic Fitting: 


Selec uta | data 1 x File Edit View Insert Tools Desktop Window Help eT 
| (| Center and scale x data JO sa §/8Qnrs|\2\/08)|sa 
Plot fits Numerical results } 
| Check to display fits on figure ar 7 | 
) spline interpolant (a Fit. [linear = | 
(1) shape-preserving interpolant Coefficients and norm of residuals 
a] linear y = pl*x + p2 
1) quadratic 
©) cubic = Coefficients: 
[) 4th degree polynomial pl = 1.5584 
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1) 6th degree polynomial 
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As shown above, the relationship between temperature and internal 
energy is: 


y = 1.62 + 2347.2 
or 
internal energy = 1.6 temperature + 2347.2 


Exercise: 
Problem: 
Using the following experimental values [footnote], plot a velocity- 


time graph and determine the equation, relating the velocity and time 
for a moving object. 


Velocity [m/s] Time [s] 


12 0 


Velocity [m/s] 
142 

912 

1122 

1972 


Experimental data. 


Time [s] 
5 

10 

15 


20 


Engineering Science by E. Hughes and C. Hughes, Longman © 1994, 


(p. 375) 


Solution: 


Data can be entered as follows: velocity=[12 142 512 1122 
1972]; time=[0 5 10 15 20]; wecannow plot the data by 
typing in plot(time, velocity) ;title('Velocity-Time 
Graph');xlabel('time');ylabel('velocity' ); atthe 
MATLAB prompt. The following plot is generated, select Tools > 
Basic Fitting, notice that we are choosing the quadratic option this 


time: 

A Basic Fitting - 1 =|6)| x BD Figure 1 el| 2% 
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As shown above, the relationship between velocity and time is: 


y = 4.827 + 27 +12 


Generating Reports with MATLAB 
Basic Publishing with MATLAB 


Fablishing with MATLAB 
Commentary on the Code 
MATLAB Code 


Results of the executed Code 


e 


www, hetemeal.com 


MATLAB includes an automatic report generator called publisher. The 
publisher publishes a script in several formats, including HTML, XML, MS 
Word and PowerPoint. The published file can contain the following: 


e Commentary on the code, 

e MATLAB code, 

e Results of the executed code, including the Command Window output 
and figures created by the code. 


The 
publish 


Function 


The most basic syntax is publish('file', 'format' ) where the m- 
file is called and executed line by line then saved to a file in specified 
format. All published files are placed in the html directory although the 
published output might be a doc file. 


Publishing with Editor 


The publisher is easily accessible from the Tool Strip: 


a MATLAB R2015a Ke = 


SS (iSearch Documentation SiR 


Publish tab on the Tool Strip. 


a Edit Configurations 


publishing.m *' Publish configuration name: publishing 


= publishing.m 
oO publishing MATLAB expression: 


\* Modify expression to add input arguments. 
|= Example: 

id a=[12 3; 45 6]7 

|* foo(a); 


| 
jpublishing 


Publish settings: User Default v | Save As... 

le Output settings ies 
Output file format html 
Output folder H:\MATLAB\2015\html 
XSL file 


= Figure settings 
Figure capture method entireGUIWindow 
Image Format default (png) 
Use new figure true 
| Max image width (pixels) Inf 
Max image height (pixels) Inf 


Create thumbnail true 
=| Code settings 
Include code true 
Evaluate code true ¥ 


Close Publish Help 


Editing publishing options, currently output format is html. 


Example: 

Write a simple script and publish it in an html file. 

Select File > New > Script to create an m-file. Once the editor is opened, 
type in the following code: 


x [0:6]; % Create a row vector 

y = 1.6*x; % Compute y as a function of x 

pee ace dl % Transpose vectors x and y 
plot(x,y),title('Graph of 
y=f(x)'),xlabel('x'),ylabel('f(x)'),grid % Plot a 
graph 


Save the script as publishing.m and select Publish > Publish. 


4 MATLAB R2015b - academic 


> FE b> c > Users»! Documer 


Command Window 


5 
3 
(2 New to MATLAB? See resources for Getting Started. 
E fe >> is x = [0:6]; % Create a row vector 
3 Bhi y = 1.6*x; % Compute y as a function of x 
aT Bie Iky') % Transpose vectors x and y 
4- plot (x,y),title('Graph of y=f(x)'),xlabel('x'),ylabel('f(x)'),grid % Plot a graph 


Publishing a script. 


An HTML file is generated as shown in the figure below: 

@ Web Browser - publishing [ea |tal| 3% | 
‘J publishing x | + | 0 (0) ~ 
enc i) TS G = 


x = [0:6];  % Create a row vector 
y = 1.6*x; % Compute y as a function of x 
[x',y'] % Transpose vectors x and y 


plot (x,y) ,title('Graph of y=£(x)'),xlabel('x'),ylabel('f(x)'),grid % Plot a graph 


ans = 


Graph of y=f(x) 


A script published in html 


The Double Percentage 
%% 
Sign 


The scripts sometimes can be very long and their readability might be 
reduced. To improve the publishing result, sections are introduced by 
adding descriptive lines to the script preceded by %%. Consider the 
following example. 


Example: 
Edit the script created in the example above to look like the code below: 


%% This file creates vectors, displays results and 
plots an x-y graph 

xX = [0:6]; % Create a row vector 

y = 1.6*x; % Compute y as a function of x 

%% Tabulated data 

ix, y¥ | % Transpose vectors x and y 

%% Graph of y=f(x) 

plot(x,y),title('Graph of 
y=f(x)'),xlabel('x'),ylabel('f(x)'),grid % Plot a 
graph 


Save the script, anew HTML file is generated as shown in the figure 
below: 


@ Web Browser - publishit = 
publishing x | + | & 0 8{0) ~ 
eel sia | Location: fitey//CUsers/ ARE Oocuments/MATLAB/2tm/publishing htm v 


Contents 
= This file creates vectors, displays results and plots an x-y graph 
= Tabulated data 
= Graph of y=f(x) 


This file creates vectors, displays results and plots an x-y graph 


x= [0:6]; § Create a row vector 
y= 1.6*"; % Compute y as a function of x 


Tabulated data 


[x',y"] % Transpose vectors x and y 


Graph of y=f(x) 


plot(x,y),title('Gzaph of y=£(x)'),xlabel('x'),ylabel('=(x)'),gzid % Plot @ graph 


Graph of y=f(x) 


10 
9 aa 
8 


— — —— 


An html file with sections 


Summary of Key Points 


1. MATLAB can generate reports containing commentary on the code, 
MATLAB code and the results of the executed code, 

2. The publisher generates a script in several formats, including HTML, 
XML, MS Word and PowerPoint. 

3. The Double Percentage %% can be used to creates hyper-linked 


sections. 


Problem Set 
Problem Set for Publishing with MATLAB 
Exercise: 


Problem: 


Write a script to plot function y = aint) for 79 < @ < 107 using 


increments of 759. Publish your m-file to html. 
Solution: 


The m-file content: % This script plots a graph of 
Graph of y=sin(x)/x clc % Clear screen x = 
pi/100:pi/100:10*pi; % Create a row vector y = 
Sim(je/x: 7% Galculate y as Tunction OF x 
plot(x,y),trtle( Graph of 
y=sin(x)/x'),xlabel('x'),ylabel('y'), grid The 
html output: 


H:/MATLAB/EngineeringComputationWithMATLAB/html/Publish0: 


% This script plots a graph of Graph of y=sin(x)/x 

cle % Clear screen 

x = pi/lOO:pi/100:10*pi; % Create a row vector 

y = sin(x)./x; % Calculate y as function of x 

plot (x,y),title('Graph of y=sin(x)/x'),xlabel('x'),ylabel('y'),grid 


Graph of y=sinGe/x 


The published html file. 


Exercise: 


Problem: 


A gas is expanded in an engine cylinder, following the law PV1.3=c. 
The initial pressure is 2550 kPa and the final pressure is 210 kPa. If the 
volume at the end of expansion is 0.75 m3, write a script to compute 
the work done by the gas and publish your solution to an html file. 
This is the same problem as this Problem you have solved before. 


Solution: 


The m-file content: cle disp('A gas is expanded in an 
engine cylinder, following the law PV‘1.3=c') 
disp('The initial pressure is 2550 kPa and the 
final pressure is 210 kPa.') disp('If the 
volume at the end of expansion is 0.75 m3,') 
disp('Compute the work done by the gas.') 
disp(' ') % Display blank line n=1.3; P_1=2550; 
% Initial pressure P_f=210; % Final pressure 
V_f=.75; % Final volume V_i=(P_f* 
(V_f4n)/P_1)4(1/n); % Initial volume 
c=P_f*V_f4n; v=V_1i:.001:V_f; % Creating a row 
vector for volume, v p=c./(v.4n); % Computing 
pressure for volume WorkDone=trapz(v,p) % 
Integrating p*dv The html output: 


\TLAB/EngineeringComputationWithMATLAB/html/PublishO2.html 


ele 

disp('A gas is expanded in an engine cylinder, following the law PV*1.3=c') 
disp('The initial pressure is 2550 kPa and the final pressure is 210 kPa.') 
disp('If the volume at the end of expansion is 0.75 m3,') 

disp('Compute the work done by the gas.') 

disp(' ') % Display blank line 

n=1.3; 

P_i=2550; % Initial pressure 

P_£=210; % Final pressure 

Wage t eg % Final volume 

V_i=(P_f£*(V_f*n)/P_i)*(1/n); %* Initial volume 

a ae bee 

v=V_i:.001:V_f; % Creating a row vector for volume, v 

p=c./(v.*n); % Computing pressure for volume 

WorkDone=trapzi(v,p) % Integrating p*dv 


A gas is expanded in an engine cylinder, following the law PV*1.3=c 
The initial pressure is 2550 kPa and the final pressure is 210 kPa. 
If the volume at the end of expansion is 0.75 m3, 

Compute the work done by the gas. 


WorkDone = 


409.0666 


Published with MATLABB 7.1t 


The published html file. 


Exercise: 


Problem: 


A force F acting on a body at a distance s from a fixed point is given 
by fF = 3s+ +. Write a script to compute the work done when the 
body moves from the position where s=1 to that where s=10 and and 
publish your solution to an html file. Include a table of contents in the 
output file. This is the same problem as this Problem you have solved 
before. 


Solution: 


The m-file content: %% Work done % This script 
computes the work done on an object clc disp('A 
force F acting on a body at a distance s froma 
fixed point is given by') disp('F=3*s+(1/(s42) ) 
where s is the distance in meters' ) 
disp('Compute the total work done in moving' ) 
disp('From the position where s=1 to that where 
s=10.') disp(' ') % Display blank line %% 
Create a row vector for distance, s 
S=1:.001:10; %% Compute Force for s F=3.*S+ 
(1./(s.42)); % Computing Force for s %% 
Integrating F*ds over 1 to 10 meters. 
WorkDone=trapz(s,F) The html output: 


\TLAB/EngineeringComputationWithMATLAB/html/PublishO3.html 


Contents 


= Work done 
« Create a row vector for distance, s 


= Compute Force for s 
= Integrating F*ds over 1 to 10 meters. 


Work done 


This script computes the work done on an object 


cle 

disp('A force F acting on a body at a distance s from a fixed point is given by') 
disp('F=3*s+(1/ (s*2)) where s is the distance in meters') 

disp('Compute the total work done in moving') 

disp('From the position where s=l1 to that where s=10.') 

disp(' ') % Display blank line 


A force F acting on a body at a distance s from a fixed point is given by 
F=3*s+(1l/(s*2)) where s is the distance in meters 

Compute the total work done in moving 

From the position where s=l1 to that where s=10. 


Create a row vector for distance, s 

s=1: .001:10; 
Compute Force fors 

F=3.*s+(1./(s.*2)); % Computing Force for s 
Integrating F*ds over 1 to 10 meters. 


VorkDone=trapz(s,F) 


WorkDone = 


149.4000 


The published html file. 


